BACKGROUND
analysis of genetic variation in human populations is critical to the understanding of the genetic basis for complex diseases.
most studied differences in dna are single-nucleotide variations at particular positions in the genome, which are called single nucleotide polymorphisms .
the positions are also called markers and the different possible values alleles.
a haplotype is a sequence of snp alleles along a region of a chromosome, and concisely represents the  genetic information in that region.
in the search for dna sequence variants that are related to common diseases , haplotype-based approaches have become a central theme  <cit> .
in diploid organisms such as humans there are two homologous  copies of each chromosome.
current practical laboratory measurement techniques produce a genotype – for m markers, a sequence of m unordered pairs of alleles.
the genotype reveals which two alleles are present at each marker, but not their respective chromosomal origin.
in order to obtain haplotypes from genotype data, this hidden phase information needs to be reconstructed.
there are two alternative approaches: if family trios are available, most of the ambiguity in the phase can be resolved analytically.
if not, population-based computational methods have to be used to estimate the haplotype pair for each genotype.
because trios are more difficult to recruit and more expensive to genotype, population-based approaches are often the only cost-effective method for large-scale studies.
consequently, the study of such techniques has received much attention recently  <cit> .
in this paper, we propose and evaluate a novel approach for population-based haplotyping based on constrained hidden markov models.
population-based haplotype reconstruction
a haplotype h is a sequence of alleles h in markers i =  <dig> ...,m. in most cases, only two alternative alleles occur at an snp marker, so we can assume that h ∈ { <dig>  1}m. a genotype g is a sequence of unordered pairs g = {hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@} of alleles in markers i =  <dig> ...,m. hence, g ∈ {{ <dig>  0}, { <dig>  1}, { <dig>  1}}m. a marker with alleles { <dig>  0} or { <dig>  1} is homozygous whereas a marker with alleles { <dig>  1} is heterozygous.
problem  <dig>
given a multiset g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ of genotypes, find for each g ∈ g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ the most likely haplotypes hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@ and hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@ which are a consistent reconstruction of g, i.e., g = {hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@} for each i =  <dig> ...,m.
if ℋ
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ denotes a mapping g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ → { <dig>  1}m × { <dig>  1}m, associating each genotype g ∈ g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ with a pair ⟨hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩ of haplotypes, the goal is to find the ℋ
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ that maximizes p. it is usually assumed that the sample g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ is in hardy-weinberg equilibrium, i.e., that p = pp for all g ∈ g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@, and that genotypes are independently sampled from the same distribution.
with such assumptions, the likelihood p of the reconstruction ℋ
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ given g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ is proportional to ∏g∈gpp
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqeqaqaamrr1ngbprwthrhayaqeguudjxwakbsthrhagq1dvbaceagae8xghalaeiikagiaemiaag2aa0baasqaaiabdeganbqaaiabigdaxaaakiabcmcapiab=lriqjabcicaoiabdigaonaadaaaleaacqwgnbwzaeaacqaiyagmaagccqggpaqkasqaaiabdeganjabgigioprthrhal1wy0l2yhvtyaexbnfgdovwbhrxajfwnagqbaiab+zq8hbqab0gaey4diunaaaa@53ae@ if the reconstruction is consistent for all g ∈ g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@, and zero otherwise.
in population-based haplotyping, a probabilistic model λ for the distribution over haplotypes is estimated from the available genotype information g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@.
the distribution estimate p is then used to find the most likely reconstruction ℋ
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ for g
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ under hardy-weinberg equilibrium.
the genetic variation in snps is mostly due to two causes: mutation and recombination.
mutations are relatively rare, they occur with a frequency of about 10- <dig>  while snps are themselves results of ancient mutations, mutations are usually ignored in statistical haplotype models due to their rarity.
recombination introduces variability by breaking up the chromosomes of the two parents and reconnecting the resulting segments to form a new and different chromosome for the offspring.
because the probability of a recombination event between two markers is lower if they are near to each other, there is a statistical correlation  between markers which decreases with increasing marker distance.
statistical approaches to haplotype modeling are based on exploiting such patterns of correlation.
methods
this section presents the proposed method for haplotype reconstruction.
we discuss the statistical model employed and present an incremental algorithm for efficiently learning the model structure from genotype data.
finally, datasets and systems used in the experimental evaluation are described.
markov models for haplotyping
we model the probability distribution on haplotypes by a left-right markov model λ with 2·m states, with a state space as shown in figure  <dig>  a haplotype  is sampled by traversing a path through the model from left to right.
the markov assumption p=∏t=1mpt
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiab=lriqjabcicaoiabdigaojabcmcapiabg2da9maaradabagae8xgha1aasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabcqggpaqkasqaaiabdsha0jabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@571c@ is motivated by the observation that linkage disequilibrium decreases with increasing marker distance.
parameters are of the form pt, the probability of sampling the new allele h at position t after observing the allele h at position t -  <dig>  note that separate  allele distributions pt are defined for every sequence position t ∈ { <dig> ...,m}, as linkage disequilibrium patterns will vary for different markers.
this also means that the allele encoding at a given marker position, i.e., which allele is represented as '0' and which as '1', does not affect the distributions that can be represented.
this model is not directly applicable in haplotype reconstruction, because in reality only genotypes are observed whereas the phase information is hidden.
the hidden phase information can be modeled by a hidden markov model λ' as shown in figure  <dig>  a path through this model corresponds to sampling a pair of haplotypes , while the corresponding genotype  is emitted.
to reflect the hardy-weinberg equilibrium assumption, constraints have to be placed on transition probabilities.
a transition in this model corresponds to independently sampling two new alleles h <dig> and h <dig> at marker t based on their respective histories h <dig> and h <dig>  therefore, the corresponding probability is actually the product of probabilities for sampling hi after hi:
pt = ptpt.
in this way, all parameters of λ' can be re-expressed as products of parameters of the model λ on haplotypes outlined above.
furthermore, λ' can be transformed into an equivalent hmm in which these constraints involving products of parameters are replaced with standard parameter tying constraints, which tie parameters in λ' to those in λ.
an advantage of this approach is that the model λ' can be trained directly from genotype data using baum-welsh algorithm  <cit> , while implicitly estimating the distribution over haplotypes encoded in λ. furthermore, the most likely reconstruction of a genotype can be directly obtained by the viterbi algorithm  <cit> .
the presented idea of embedding a model on haplotypes into a model on genotypes in which the genotype phase is the hidden state information, and learning this model using em, is related to the approaches used in the hit  <cit>  and fastphase  <cit>  systems.
in hit, haplotypes are modeled as recombinations of a set of founder haplotypes, and an instance of the em algorithm is derived to directly estimate the founders from genotype observations.
in fastphase, haplotypes are modeled using local clusters, and cluster membership of a haplotype is determined by a hidden markov model.
again, an instance of the em algorithm for estimating the clusters directly from genotype data can be derived.
higher-order models and sparse distributions
the main limitation of the model presented so far is that it only takes into account dependencies between adjacent markers.
expressivity can be increased by using a markov model of order k >  <dig> for the underlying haplotype distribution  <cit> :
p=∏t=1mpt,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiab=lriqjabcicaoiabdigaojabcmcapiabg2da9maarahabagae8xgha1aasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabcqggpaqkcqggsaalasqaaiabdsha0jabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@5cd9@
where h is a shorthand for h...h. unfortunately, the number of parameters in such a model increases exponentially with the history length k. fortunately, observations on real-world data  show that only few conserved haplotype fragments from the set of 2k possible binary strings of length k actually occur.
this can be exploited by modeling sparse distributions, where fragment probabilities which are estimated to be very low are set to zero.
more precisely, let p = pt and define for some small ε >  <dig> a regularized distribution
p^t={0if p≤ε;1ifp>1−ε;potherwise.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaajawaasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggpaqkcqgh9aqpdagabaqaauaabaaadiaaaeaacqaiwaamaeaaieaacqgfpbqacqgfmbgzcqqggaaicqwgwbaccqghkjyoiigacqqf1oqzcqgg7awoaeaacqaixaqmaeaacqgfpbqacqgfmbgztcvaufebsjuyzl2yd9gzlbvynv2caexbulwblnhiov2dgi1btfmbagqbaiaa8bcacqwgwbaccqgh+agpcqaixaqmcqghsislcqqf1oqzcqgg7awoaeaacqwgwbacaeaacqgfvbwbcqgf0badcqgfobaacqgflbqzcqgfybgccqgf3bwdcqgfpbqacqgfzbwccqgflbqzcqgguaglaaaacagl7baaaaa@80ab@
if the underlying distribution is sufficiently sparse, p^
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaajaaaaa@376d@ can be represented using a relative small number of parameters.
the corresponding sparse markov model structure  will reflect the pattern of conserved haplotype fragments present in the population.
how such a sparse model structure can be learned without ever constructing the prohibitively complex distribution p will be discussed in the next section.
spamm: a level-wise learning algorithm
algorithm  <dig> the level-wise spamm learning algorithm.
initialize k := 1
λ <dig> := initial-model()
λ <dig> := em-training
repeat
k := k + 1
λk := extend-and-regularize
λk := em-training
until k = kmax
to construct the sparse order-k hidden markov model, we propose a learning algorithm – called spamm for sparse markov modeling – that iteratively refines hidden markov models of increasing order .
more specifically, the idea of spamm is to identify conserved fragments using a level-wise search, i.e., by extending short fragments  to longer ones , and is inspired by the well-known apriori data mining algorithm  <cit> .
the algorithm starts with a first-order markov model λ <dig> on haplotypes where initial transition probabilities are set to p˙
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaacaaaaa@3766@t =  <dig>  for all t ∈ { <dig> ...,m}, h, h ∈ { <dig>  1}.
this model can be embedded into a hidden markov model λ′1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ on genotypes as explained above, and λ′1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ can be trained from the available genotype data using the standard em algorithm.
as parameters in λ′1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ are tied to those in λ <dig>  this yields new estimates for the parameters pt in λk.
this training procedure is summarized in the function em-training.
the function extend-and-regularize takes as input a model of order k -  <dig> and returns a model λk of order k. in λk, initial transition probabilities are set to
p˙t≤ε;1if pt>1−ε; <dig> otherwise,
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaacawaasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabdawgaawcbagaem4aasmaey4kasiaegymaedabeaakiabg2da9maaceaabaqbaeaaamgaaaqaaiabicdawaqaagqaaiab9lgapjab9zgamjaaykw7cqwfzecudawgaawcbagaemidaqhabeaakiabcicaoiabdigaojabcufabjabdsha0jabc2fadjabcyha8jabdigaojabcufabjabdsha0jabgkhitiabdugarjabgucariabigdaxiabcycasiabdsha0jabgkhitiabigdaxiabc2fadjabcycasiab+t7asnaabaaaleaacqwgrbwaaeqaaogaeiykakiaeyizimqae4xtdumaei4oasdabagaegymaedabagae0xaakmae0nzaymaagpavlab=lriqnaabaaaleaacqwg0badaeqaaogaeiikagiaemiaagmaei4waslaemidaqnaeiyxa0laeiifawnaemiaagmaei4waslaemidaqnaeyoei0iaem4aasmaey4kasiaegymaejaeiilawiaemidaqnaeyoei0iaegymaejaeiyxa0laeiilawiae43udw2aasbaasqaaiabdugarbqabagccqggpaqkcqgh+agpcqaixaqmcqghsislcqgf1oqzcqgg7awoaeaacqaiwaamcqgguaglcqai1aqnaeaacqqfvbwbcqqf0badcqqfobaacqqflbqzcqqfybgccqqf3bwdcqqfpbqacqqfzbwccqqflbqzcqggsaalaaaacagl7baaaaa@b291@
i.e., transitions are removed if the probability of the transition conditioned on a shorter history is smaller than ε. this procedure of iteratively training, extending and regularizing markov models of increasing order is repeated up to a maximum order kmax.
for a given genotype g, a reconstructed haplotype pair ⟨hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩k can be obtained from every model λk.
at the same time, the viterbi algorithm computes p, an estimate of the confidence of the reconstruction.
in spamm, the reconstruction ⟨hg1
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩k* with the highest confidence is returned as the final solution:
k∗=arg⁡max⁡k∈{ <dig> ...,kmax}p.
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgrbwacqghxiikcqgh9aqpdawfqaqaaigbcggahjabckhayjabceganjgbc2gatjabcggahjabciha4bwcbagaem4aasmaeyici4saei4easnaegymaejaeiilawiaeiola4iaeiola4iaeiola4iaeiilawiaem4aas2aasbaawqaagqaciab=1gatjab=fgahjab=hha4bqabawccqgg9bqfaeqaamrr1ngbprwthrhayaqeguudjxwakbsthrhagq1dvbaceagccqgfzecucqggoaakcqghpms4cqwgobaadaqhaawcbagaem4zacgabagaegymaedaaogaeiilawiaemiaag2aa0baasqaaiabdeganbqaaiabikdayaaakiabgqyixpaabaaaleaacqwgrbwaaeqaaogaeiifawnaem4zacmaeiilawcccigae03udw2aasbaasqaaiabdugarbqabagccqggpaqkcqgguaglaaa@6b66@
the idea of using frequent fragments to build markov models for haplotypes has also been used in the haplorec method  <cit> .
in haplorec, a set of fragments  that are frequent according to the current model is kept, and updated after each iteration of the em algorithm.
experimental methodology and evaluation
the proposed method was implemented in the spamm haplotyping system  <cit> .
we compared its accuracy and computational performance to several other state-of-the art haplotype reconstruction systems: phase version  <dig> .
<dig>  <cit> , fastphase version  <dig>   <cit> , gerbil as included in gevalt version  <dig>   <cit> , hit  <cit>  and haplorec  version  <dig>   <cit> .
all methods were run using their default parameters.
the fastphase system, which also employs em for learning a probabilistic model, uses a strategy of averaging results over several random restarts of em from different initial parameter values.
this reduces the variance component of the reconstruction error and alleviates the problem of local minima in em search.
as this is a general technique applicable also to our method, we list results for fastphase with averaging  and without averaging .
the methods were compared using publicly available real-world datasets, and larger datasets simulated with the hudson coalescence simulator  <cit> .
as real-world data, we used a collection of datasets from the yoruba population in ibadan, nigeria  <cit> , and the well-known dataset of daly et al  <cit> , which contains data from a european-derived population.
for these datasets, family trios are available, and thus true haplotypes can be inferred analytically.
non-transmitted parental chromosomes of each trio were combined to form additional artificial haplotype pairs.
markers with minor allele frequency of less than 5% and genotypes with more than 15% missing values were removed.
note that if all trio members are heterozygous, the haplotype of the child can not be inferred.
in this case, the genotype at this marker position is observed but the marker is ignored when computing the accuracy of the method.
for the yoruba population, information on  <dig>  million snps spread over the whole genome is available.
we sampled  <dig> sets of  <dig> markers each from distinct regions on chromosome  <dig> , and from these smaller datasets by taking only the first  <dig>  or  <dig>  markers for every individual.
there are  <dig> individuals in the dataset after preprocessing, with an average fraction of missing values of  <dig> %.
for the daly dataset, there is information on  <dig> markers and  <dig> individuals available after data preprocessing, and the average fraction of missing values is 8%.
although results on a single dataset are not very meaningful, the daly dataset was included because it has been used frequently in the literature.
the number of genotyped individuals in these real-world datasets is rather small.
for most disease association studies, sample sizes of at least several hundred individuals are needed  <cit> , and we are ultimately interested in haplotyping such larger datasets.
unfortunately, we are not aware of any publicly available real-world datasets of this size, so we have to resort to simulated data.
we used the well-known hudson coalescence simulator  <cit>  to generate  <dig> artificial datasets, each containing  <dig> individuals .
the simulator uses the standard wright-fisher neutral model of genetic variation with recombination.
a chromosomal region of  <dig> kb was simulated.
the probability of mutation in each base pair was set to 10- <dig> per generation, and the probability of cross-over between adjacent base pairs was set to 10- <dig>  these values result in a mutation probability for the entire chromosomal region of μ =  <dig>  and cross-over probability of ρ =  <dig> .
the diploid population size, n <dig>  was set to the standard  <dig>  yielding mutation parameter θ = 4n0μ =  <dig>  and the recombination parameter r =  <dig>  for each data set, a sample of  <dig> chromosomes was generated, and these were paired to form  <dig> genotypes.
on average, one simulation produced approximately  <dig> segregating sites.
for each data set,  <dig> markers were chosen from the segregating sites with minor allele frequency of at least 5%, such that marker spacing was as uniform as possible.
the resulting average marker spacing was  <dig>  kb.
to come as close to the characteristics of real-world data as possible, some alleles were masked  after simulation.
more specifically, the missing allele pattern found in the yoruba datasets was superimposed onto the simulated data, shortening patterns to the size of the target marker map and repeating them as needed for additional individuals.
the accuracy of the reconstructed haplotypes produced by the different methods was measured by normalized switch error.
the switch error of a reconstruction is the minimum number of recombinations needed to transform the reconstructed haplotype pair into the true haplotype pair.
to normalize, switch errors are summed over all individuals in the dataset and divided by the total number of switch errors that could have been made.
RESULTS
performance of the systems on larger datasets with up to  <dig> individuals was evaluated on the  <dig> simulated hudson datasets.
as for the real-world data, the most accurate methods were phase, fastphase, spamm and haplorec.
figure  <dig> shows the normalized switch error of these four methods as a function of the number of individuals .
phase was the most accurate method also in this setting, but the relative accuracy of the other three systems depended on the number of individuals in the datasets.
while for relatively small numbers of individuals  fastphase outperforms spamm and haplorec, this is reversed for  <dig> or more individuals.
a problem closely related to haplotype reconstruction is that of genotype imputation.
here, the task is to infer the most likely genotype values  at marker positions where genotype information is missing, based on the observed genotype information.
with the exception of haplorec, all haplotyping systems included in this study can also impute missing genotypes.
to test imputation accuracy, between 10% and 40% of all markers were masked randomly, and then the marker values inferred by the systems were compared to the known true marker values.
table  <dig> shows the accuracy of inferred genotypes for different fractions of masked data on the yoruba- <dig> datasets and table  <dig> on the simulated hudson datasets with  <dig> individuals per dataset.
phase was too slow to run in this task as its runtime increases significantly in the presence of many missing markers.
evidence from the literature  <cit>  suggests that for this task, fastphase outperforms phase and is indeed the best method available.
in our experiments, on yoruba- <dig> fastphase is most accurate, spamm is slightly less accurate than fastphase, but more accurate than any other method .
on the larger hudson datasets, spamm is significantly more accurate than any other method.
our experimental results confirm phase as the most accurate but also computationally most expensive haplotype reconstruction system  <cit> .
if more computational efficiency is required, fastphase yields the most accurate reconstructions on small datasets, and spamm is preferable for larger datasets.
spamm also infers missing genotype values with high accuracy.
for small datasets, it is second only to fastphase; for large datasets, it is substantially more accurate than any other method in our experiments.
the presented method is quite basic: it does not use fine-tuned priors for em, multiple em restarts or averaging techniques  <cit> , or cross-validates model parameters  <cit> .
moreover, most statistical models employed in haplotyping are specifically tailored to this problem, and reflect certain assumptions about haplotype structure.
for example, the hit method assumes that there is a limited number of founder haplotypes for a population, and gerbil assumes block-like haplotype patterns.
these systems are only effective if the underlying assumptions are valid.
hit, for instance, was less accurate than phase in our study, but has been shown to be competitive with phase on population samples from finland  <cit> , a population isolate for which the assumption of a small number of founders is particularly realistic  <cit> .
similarly, performance of gerbil will suffer if haplotypes do not exhibit a block-like structure.
in contrast, the sparse higher-order markov chains used in spamm are a general sequence modeling technique.
detailed assumptions about the haplotype structure are replaced by the structure-learning component of the algorithm.
the resulting model is rather flexible, and subsumes block-like or mosaic-like haplotype structures .
in fact, the proposed approach is not limited to haplotype analysis, and an interesting direction for future work is to apply it also to other sequence modeling tasks.
CONCLUSIONS
we proposed a simple haplotype reconstruction method that is based on iterative refinement and regularization of constrained hidden markov models .
the method was compared against several other state-of-the-art haplotyping systems on real-world genotype datasets with 60– <dig> individuals and larger simulated datasets with up to  <dig> individuals.
in the experimental study, phase was the most accurate, but also computationally most demanding haplotype reconstruction system.
fastphase and spamm are slightly less accurate but much faster, and scale well to long marker maps.
the relative performance of these two systems depends on the number of samples available: while fastphase is slightly more accurate for small datasets, spamm is superior for datasets with several hundred genotype samples.
as large datasets are ultimately needed for successful disease association studies, the presented method is a promising alternative to existing approaches.
authors' contributions
tm, nl and hm developed the haplotyping method.
nl implemented the method and carried out the experiments.
le contributed data to the experimental evaluation.
hm and ht coordinated the research.
all authors contributed to the preparation of the manuscript.

 analysis of genetic variation in human populations is critical to the understanding of the genetic basis for complex diseases. most studied differences in dna are single-nucleotide variations at particular positions in the genome, which are called single nucleotide polymorphisms. the positions are also called markers and the different possible values alleles. a haplotype is a sequence of snp alleles along a region of a chromosome, and concisely represents the genetic information in that region. in the search for dna sequence variants that are related to common diseases, haplotype-based approaches have become a central theme. in diploid organisms such as humans there are two homologous copies of each chromosome. current practical laboratory measurement techniques produce a genotype – for m markers, a sequence of m unordered pairs of alleles. the genotype reveals which two alleles are present at each marker, but not their respective chromosomal origin. in order to obtain haplotypes from genotype data, this hidden phase information needs to be reconstructed. there are two alternative approaches: if family trios are available, most of the ambiguity in the phase can be resolved analytically. if not, population-based computational methods have to be used to estimate the haplotype pair for each genotype. because trios are more difficult to recruit and more expensive to genotype, population-based approaches are often the only cost-effective method for large-scale studies. consequently, the study of such techniques has received much attention recently. in this paper, we propose and evaluate a novel approach for population-based haplotyping based on constrained hidden markov models. population-based haplotype reconstruction a haplotype h is a sequence of alleles h in markers i =...,m. in most cases, only two alternative alleles occur at an snp marker, so we can assume that h ∈ { 1}m. a genotype g is a sequence of unordered pairs g = {hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@} of alleles in markers i =...,m. hence, g ∈ {{ 0}, { 1}, { 1}}m. a marker with alleles { 0} or { 1} is homozygous whereas a marker with alleles { 1} is heterozygous. problem given a multiset g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ of genotypes, find for each g ∈ g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ the most likely haplotypes hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@ and hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@ which are a consistent reconstruction of g, i.e., g = {hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@} for each i =...,m. if ℋ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ denotes a mapping g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ → { 1}m × { 1}m, associating each genotype g ∈ g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ with a pair ⟨hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩ of haplotypes, the goal is to find the ℋ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ that maximizes p. it is usually assumed that the sample g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ is in hardy-weinberg equilibrium, i.e., that p = pp for all g ∈ g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@, and that genotypes are independently sampled from the same distribution. with such assumptions, the likelihood p of the reconstruction ℋ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ given g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ is proportional to ∏g∈gpp mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqeqaqaamrr1ngbprwthrhayaqeguudjxwakbsthrhagq1dvbaceagae8xghalaeiikagiaemiaag2aa0baasqaaiabdeganbqaaiabigdaxaaakiabcmcapiab=lriqjabcicaoiabdigaonaadaaaleaacqwgnbwzaeaacqaiyagmaagccqggpaqkasqaaiabdeganjabgigioprthrhal1wy0l2yhvtyaexbnfgdovwbhrxajfwnagqbaiab+zq8hbqab0gaey4diunaaaa@53ae@ if the reconstruction is consistent for all g ∈ g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@, and zero otherwise. in population-based haplotyping, a probabilistic model λ for the distribution over haplotypes is estimated from the available genotype information g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@. the distribution estimate p is then used to find the most likely reconstruction ℋ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwflecsaaa@3763@ for g mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbamrthrhal1wy0l2yhvtyaehbnfgdovwbhrxajfwnaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabawaaegaeaaakeaaimaacqwfge=raaa@382d@ under hardy-weinberg equilibrium. the genetic variation in snps is mostly due to two causes: mutation and recombination. mutations are relatively rare, they occur with a frequency of about 10- while snps are themselves results of ancient mutations, mutations are usually ignored in statistical haplotype models due to their rarity. recombination introduces variability by breaking up the chromosomes of the two parents and reconnecting the resulting segments to form a new and different chromosome for the offspring. because the probability of a recombination event between two markers is lower if they are near to each other, there is a statistical correlation between markers which decreases with increasing marker distance. statistical approaches to haplotype modeling are based on exploiting such patterns of correlation. methods this section presents the proposed method for haplotype reconstruction. we discuss the statistical model employed and present an incremental algorithm for efficiently learning the model structure from genotype data. finally, datasets and systems used in the experimental evaluation are described. markov models for haplotyping we model the probability distribution on haplotypes by a left-right markov model λ with 2·m states, with a state space as shown in figure a haplotype is sampled by traversing a path through the model from left to right. the markov assumption p=∏t=1mpt mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiab=lriqjabcicaoiabdigaojabcmcapiabg2da9maaradabagae8xgha1aasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabcqggpaqkasqaaiabdsha0jabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@571c@ is motivated by the observation that linkage disequilibrium decreases with increasing marker distance. parameters are of the form pt, the probability of sampling the new allele h at position t after observing the allele h at position t - note that separate allele distributions pt are defined for every sequence position t ∈ {...,m}, as linkage disequilibrium patterns will vary for different markers. this also means that the allele encoding at a given marker position, i.e., which allele is represented as'0' and which as'1', does not affect the distributions that can be represented. this model is not directly applicable in haplotype reconstruction, because in reality only genotypes are observed whereas the phase information is hidden. the hidden phase information can be modeled by a hidden markov model λ' as shown in figure a path through this model corresponds to sampling a pair of haplotypes, while the corresponding genotype is emitted. to reflect the hardy-weinberg equilibrium assumption, constraints have to be placed on transition probabilities. a transition in this model corresponds to independently sampling two new alleles h and h at marker t based on their respective histories h and h therefore, the corresponding probability is actually the product of probabilities for sampling hi after hi: pt = ptpt. in this way, all parameters of λ' can be re-expressed as products of parameters of the model λ on haplotypes outlined above. furthermore, λ' can be transformed into an equivalent hmm in which these constraints involving products of parameters are replaced with standard parameter tying constraints, which tie parameters in λ' to those in λ. an advantage of this approach is that the model λ' can be trained directly from genotype data using baum-welsh algorithm, while implicitly estimating the distribution over haplotypes encoded in λ. furthermore, the most likely reconstruction of a genotype can be directly obtained by the viterbi algorithm. the presented idea of embedding a model on haplotypes into a model on genotypes in which the genotype phase is the hidden state information, and learning this model using em, is related to the approaches used in the hit and fastphase systems. in hit, haplotypes are modeled as recombinations of a set of founder haplotypes, and an instance of the em algorithm is derived to directly estimate the founders from genotype observations. in fastphase, haplotypes are modeled using local clusters, and cluster membership of a haplotype is determined by a hidden markov model. again, an instance of the em algorithm for estimating the clusters directly from genotype data can be derived. higher-order models and sparse distributions the main limitation of the model presented so far is that it only takes into account dependencies between adjacent markers. expressivity can be increased by using a markov model of order k > for the underlying haplotype distribution: p=∏t=1mpt, mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiab=lriqjabcicaoiabdigaojabcmcapiabg2da9maarahabagae8xgha1aasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabcqggpaqkcqggsaalasqaaiabdsha0jabg2da9iabigdaxaqaaiabd2gatbqdcqghpis1aaaa@5cd9@ where h is a shorthand for h...h. unfortunately, the number of parameters in such a model increases exponentially with the history length k. fortunately, observations on real-world data show that only few conserved haplotype fragments from the set of 2k possible binary strings of length k actually occur. this can be exploited by modeling sparse distributions, where fragment probabilities which are estimated to be very low are set to zero. more precisely, let p = pt and define for some small ε > a regularized distribution p^t={0if p≤ε;1ifp>1−ε;potherwise. mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaajawaasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggpaqkcqgh9aqpdagabaqaauaabaaadiaaaeaacqaiwaamaeaaieaacqgfpbqacqgfmbgzcqqggaaicqwgwbaccqghkjyoiigacqqf1oqzcqgg7awoaeaacqaixaqmaeaacqgfpbqacqgfmbgztcvaufebsjuyzl2yd9gzlbvynv2caexbulwblnhiov2dgi1btfmbagqbaiaa8bcacqwgwbaccqgh+agpcqaixaqmcqghsislcqqf1oqzcqgg7awoaeaacqwgwbacaeaacqgfvbwbcqgf0badcqgfobaacqgflbqzcqgfybgccqgf3bwdcqgfpbqacqgfzbwccqgflbqzcqgguaglaaaacagl7baaaaa@80ab@ if the underlying distribution is sufficiently sparse, p^ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaajaaaaa@376d@ can be represented using a relative small number of parameters. the corresponding sparse markov model structure will reflect the pattern of conserved haplotype fragments present in the population. how such a sparse model structure can be learned without ever constructing the prohibitively complex distribution p will be discussed in the next section. spamm: a level-wise learning algorithm algorithm the level-wise spamm learning algorithm.    initialize k:= 1    λ:= initial-model()    λ:= em-training    repeat       k:= k + 1       λk:= extend-and-regularize       λk:= em-training    until k = kmax to construct the sparse order-k hidden markov model, we propose a learning algorithm – called spamm for sparse markov modeling – that iteratively refines hidden markov models of increasing order. more specifically, the idea of spamm is to identify conserved fragments using a level-wise search, i.e., by extending short fragments to longer ones, and is inspired by the well-known apriori data mining algorithm. the algorithm starts with a first-order markov model λ on haplotypes where initial transition probabilities are set to p˙ mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaacaaaaa@3766@t = for all t ∈ {...,m}, h, h ∈ { 1}. this model can be embedded into a hidden markov model λ′1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ on genotypes as explained above, and λ′1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ can be trained from the available genotype data using the standard em algorithm. as parameters in λ′1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacuwf7oabgaqbamaabaaaleaacqaixaqmaeqaaaaa@2f8f@ are tied to those in λ this yields new estimates for the parameters pt in λk. this training procedure is summarized in the function em-training. the function extend-and-regularize takes as input a model of order k - and returns a model λk of order k. in λk, initial transition probabilities are set to p˙t≤ε;1if pt>1−ε; otherwise, mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaatuudjxwak1uy0hmmaehbfv3yslgzg0uy0hgiud3bagabaiqb=lriqzaacawaasbaasqaaiabdsha0bqabagccqggoaakcqwgobaacqggbbwwcqwg0badcqggdbqxcqgg8bafcqwgobaacqggbbwwcqwg0badcqghsislcqwgrbwacqggsaalcqwg0badcqghsislcqaixaqmcqggdbqxcqggsaaliigacqgf7oabdawgaawcbagaem4aasmaey4kasiaegymaedabeaakiabg2da9maaceaabaqbaeaaamgaaaqaaiabicdawaqaagqaaiab9lgapjab9zgamjaaykw7cqwfzecudawgaawcbagaemidaqhabeaakiabcicaoiabdigaojabcufabjabdsha0jabc2fadjabcyha8jabdigaojabcufabjabdsha0jabgkhitiabdugarjabgucariabigdaxiabcycasiabdsha0jabgkhitiabigdaxiabc2fadjabcycasiab+t7asnaabaaaleaacqwgrbwaaeqaaogaeiykakiaeyizimqae4xtdumaei4oasdabagaegymaedabagae0xaakmae0nzaymaagpavlab=lriqnaabaaaleaacqwg0badaeqaaogaeiikagiaemiaagmaei4waslaemidaqnaeiyxa0laeiifawnaemiaagmaei4waslaemidaqnaeyoei0iaem4aasmaey4kasiaegymaejaeiilawiaemidaqnaeyoei0iaegymaejaeiyxa0laeiilawiae43udw2aasbaasqaaiabdugarbqabagccqggpaqkcqgh+agpcqaixaqmcqghsislcqgf1oqzcqgg7awoaeaacqaiwaamcqgguaglcqai1aqnaeaacqqfvbwbcqqf0badcqqfobaacqqflbqzcqqfybgccqqf3bwdcqqfpbqacqqfzbwccqqflbqzcqggsaalaaaacagl7baaaaa@b291@ i.e., transitions are removed if the probability of the transition conditioned on a shorter history is smaller than ε. this procedure of iteratively training, extending and regularizing markov models of increasing order is repeated up to a maximum order kmax. for a given genotype g, a reconstructed haplotype pair ⟨hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩k can be obtained from every model λk. at the same time, the viterbi algorithm computes p, an estimate of the confidence of the reconstruction. in spamm, the reconstruction ⟨hg1 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegymaedaaaaa@3079@, hg2 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgobaadaqhaawcbagaem4zacgabagaegomaidaaaaa@307b@⟩k* with the highest confidence is returned as the final solution: k∗=arg⁡max⁡k∈{...,kmax}p. mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgrbwacqghxiikcqgh9aqpdawfqaqaaigbcggahjabckhayjabceganjgbc2gatjabcggahjabciha4bwcbagaem4aasmaeyici4saei4easnaegymaejaeiilawiaeiola4iaeiola4iaeiola4iaeiilawiaem4aas2aasbaawqaagqaciab=1gatjab=fgahjab=hha4bqabawccqgg9bqfaeqaamrr1ngbprwthrhayaqeguudjxwakbsthrhagq1dvbaceagccqgfzecucqggoaakcqghpms4cqwgobaadaqhaawcbagaem4zacgabagaegymaedaaogaeiilawiaemiaag2aa0baasqaaiabdeganbqaaiabikdayaaakiabgqyixpaabaaaleaacqwgrbwaaeqaaogaeiifawnaem4zacmaeiilawcccigae03udw2aasbaasqaaiabdugarbqabagccqggpaqkcqgguaglaaa@6b66@ the idea of using frequent fragments to build markov models for haplotypes has also been used in the haplorec method. in haplorec, a set of fragments that are frequent according to the current model is kept, and updated after each iteration of the em algorithm. experimental methodology and evaluation the proposed method was implemented in the spamm haplotyping system. we compared its accuracy and computational performance to several other state-of-the art haplotype reconstruction systems: phase version., fastphase version, gerbil as included in gevalt version, hit and haplorec version. all methods were run using their default parameters. the fastphase system, which also employs em for learning a probabilistic model, uses a strategy of averaging results over several random restarts of em from different initial parameter values. this reduces the variance component of the reconstruction error and alleviates the problem of local minima in em search. as this is a general technique applicable also to our method, we list results for fastphase with averaging and without averaging. the methods were compared using publicly available real-world datasets, and larger datasets simulated with the hudson coalescence simulator. as real-world data, we used a collection of datasets from the yoruba population in ibadan, nigeria, and the well-known dataset of daly et al, which contains data from a european-derived population. for these datasets, family trios are available, and thus true haplotypes can be inferred analytically. non-transmitted parental chromosomes of each trio were combined to form additional artificial haplotype pairs. markers with minor allele frequency of less than 5% and genotypes with more than 15% missing values were removed. note that if all trio members are heterozygous, the haplotype of the child can not be inferred. in this case, the genotype at this marker position is observed but the marker is ignored when computing the accuracy of the method. for the yoruba population, information on million snps spread over the whole genome is available. we sampled sets of markers each from distinct regions on chromosome, and from these smaller datasets by taking only the first or markers for every individual. there are individuals in the dataset after preprocessing, with an average fraction of missing values of %. for the daly dataset, there is information on markers and individuals available after data preprocessing, and the average fraction of missing values is 8%. although results on a single dataset are not very meaningful, the daly dataset was included because it has been used frequently in the literature. the number of genotyped individuals in these real-world datasets is rather small. for most disease association studies, sample sizes of at least several hundred individuals are needed, and we are ultimately interested in haplotyping such larger datasets. unfortunately, we are not aware of any publicly available real-world datasets of this size, so we have to resort to simulated data. we used the well-known hudson coalescence simulator to generate artificial datasets, each containing individuals. the simulator uses the standard wright-fisher neutral model of genetic variation with recombination. a chromosomal region of kb was simulated. the probability of mutation in each base pair was set to 10- per generation, and the probability of cross-over between adjacent base pairs was set to 10- these values result in a mutation probability for the entire chromosomal region of μ = and cross-over probability of ρ =. the diploid population size, n was set to the standard yielding mutation parameter θ = 4n0μ = and the recombination parameter r = for each data set, a sample of chromosomes was generated, and these were paired to form genotypes. on average, one simulation produced approximately segregating sites. for each data set, markers were chosen from the segregating sites with minor allele frequency of at least 5%, such that marker spacing was as uniform as possible. the resulting average marker spacing was kb. to come as close to the characteristics of real-world data as possible, some alleles were masked after simulation. more specifically, the missing allele pattern found in the yoruba datasets was superimposed onto the simulated data, shortening patterns to the size of the target marker map and repeating them as needed for additional individuals. the accuracy of the reconstructed haplotypes produced by the different methods was measured by normalized switch error. the switch error of a reconstruction is the minimum number of recombinations needed to transform the reconstructed haplotype pair into the true haplotype pair. to normalize, switch errors are summed over all individuals in the dataset and divided by the total number of switch errors that could have been made.  performance of the systems on larger datasets with up to individuals was evaluated on the simulated hudson datasets. as for the real-world data, the most accurate methods were phase, fastphase, spamm and haplorec. figure shows the normalized switch error of these four methods as a function of the number of individuals. phase was the most accurate method also in this setting, but the relative accuracy of the other three systems depended on the number of individuals in the datasets. while for relatively small numbers of individuals fastphase outperforms spamm and haplorec, this is reversed for or more individuals. a problem closely related to haplotype reconstruction is that of genotype imputation. here, the task is to infer the most likely genotype values at marker positions where genotype information is missing, based on the observed genotype information. with the exception of haplorec, all haplotyping systems included in this study can also impute missing genotypes. to test imputation accuracy, between 10% and 40% of all markers were masked randomly, and then the marker values inferred by the systems were compared to the known true marker values. table shows the accuracy of inferred genotypes for different fractions of masked data on the yoruba- datasets and table on the simulated hudson datasets with individuals per dataset. phase was too slow to run in this task as its runtime increases significantly in the presence of many missing markers. evidence from the literature suggests that for this task, fastphase outperforms phase and is indeed the best method available. in our experiments, on yoruba- fastphase is most accurate, spamm is slightly less accurate than fastphase, but more accurate than any other method. on the larger hudson datasets, spamm is significantly more accurate than any other method. our experimental results confirm phase as the most accurate but also computationally most expensive haplotype reconstruction system. if more computational efficiency is required, fastphase yields the most accurate reconstructions on small datasets, and spamm is preferable for larger datasets. spamm also infers missing genotype values with high accuracy. for small datasets, it is second only to fastphase; for large datasets, it is substantially more accurate than any other method in our experiments. the presented method is quite basic: it does not use fine-tuned priors for em, multiple em restarts or averaging techniques, or cross-validates model parameters. moreover, most statistical models employed in haplotyping are specifically tailored to this problem, and reflect certain assumptions about haplotype structure. for example, the hit method assumes that there is a limited number of founder haplotypes for a population, and gerbil assumes block-like haplotype patterns. these systems are only effective if the underlying assumptions are valid. hit, for instance, was less accurate than phase in our study, but has been shown to be competitive with phase on population samples from finland, a population isolate for which the assumption of a small number of founders is particularly realistic. similarly, performance of gerbil will suffer if haplotypes do not exhibit a block-like structure. in contrast, the sparse higher-order markov chains used in spamm are a general sequence modeling technique. detailed assumptions about the haplotype structure are replaced by the structure-learning component of the algorithm. the resulting model is rather flexible, and subsumes block-like or mosaic-like haplotype structures. in fact, the proposed approach is not limited to haplotype analysis, and an interesting direction for future work is to apply it also to other sequence modeling tasks.  we proposed a simple haplotype reconstruction method that is based on iterative refinement and regularization of constrained hidden markov models. the method was compared against several other state-of-the-art haplotyping systems on real-world genotype datasets with 60– individuals and larger simulated datasets with up to individuals. in the experimental study, phase was the most accurate, but also computationally most demanding haplotype reconstruction system. fastphase and spamm are slightly less accurate but much faster, and scale well to long marker maps. the relative performance of these two systems depends on the number of samples available: while fastphase is slightly more accurate for small datasets, spamm is superior for datasets with several hundred genotype samples. as large datasets are ultimately needed for successful disease association studies, the presented method is a promising alternative to existing approaches. authors' contributions tm, nl and hm developed the haplotyping method. nl implemented the method and carried out the experiments. le contributed data to the experimental evaluation. hm and ht coordinated the research. all authors contributed to the preparation of the manuscript. 
haplotype reconstruction is the problem of resolving the hidden phase information in genotype data obtained from laboratory measurements. we propose a novel approach for haplotype reconstruction based on constrained hidden markov models. relatively simple probabilistic approaches for haplotype reconstruction based on structured hidden markov models are competitive with more complex, well-established techniques in this field. results show that it is competitive with other recently proposed methods in terms of reconstruction accuracy, while offering a particularly good trade-off between computational costs and quality of results for large datasets. models are constructed by incrementally refining and regularizing the structure of a simple generative model for genotype data under hardy-weinberg equilibrium. solving this problem is an important intermediate step in gene association studies, which seek to uncover the genetic basis of complex diseases.
 haplotype reconstruction is the problem of resolving the hidden phase information in genotype data obtained from laboratory measurements. solving this problem is an important intermediate step in gene association studies, which seek to uncover the genetic basis of complex diseases. we propose a novel approach for haplotype reconstruction based on constrained hidden markov models. models are constructed by incrementally refining and regularizing the structure of a simple generative model for genotype data under hardy-weinberg equilibrium.  the proposed method is evaluated on real-world and simulated population data. results show that it is competitive with other recently proposed methods in terms of reconstruction accuracy, while offering a particularly good trade-off between computational costs and quality of results for large datasets.  relatively simple probabilistic approaches for haplotype reconstruction based on structured hidden markov models are competitive with more complex, well-established techniques in this field. 17– june probabilistic modeling and machine learning in structural and systems biology tuusula, finland 