a single dna molecule is a long chain of nucleotides ( base pairs ) . there are four such nucleotides which are represented by the set of symbols { a , t , g , c}. it is generally accepted that genomes of two humans are almost 99% identical at dna level . however , at certain specific sites , variation is observed across the human population which is commonly known as  single nucleotide polymorphism  and abbreviated as  snp  . the nucleotide involved in a snp site is called allele . if a snp site can have only two nucleotides , it is called biallelic . if it can have more than two alleles it is called a multiallelic snp . from now on , we will consider the simplest case where only bi - allelic snps occur in a specific pair of dna . the single nucleotide polymorphism ( snp ) is believed to be the most widespread form of genetic variation . the sequence of all snps in a given chromosome is called haplotype . haplotyping an individual deals with determining a pair of haplotypes , one for each copy of a given chromosome . this pair of haplotypes completely define the snp fingerprints of an individual for a specific pair of chromosomes . given the two sequences of bases , haplotyping is straight forward and just needs to iterate through both the sequences and remove all the common alleles from them . but haplotyping becomes difficult when we want to construct haplotypes from sequencing data for higher reliability . sequencing data for a genome does not contain the complete sequences of bases for a specific chromosome , rather it provides a set of fragments of arbitrary length for the whole genome . therefore , the actual problem of haplotyping is to find two haplotypes from the set of overlapping fragments of both the chromosomes , where fragments might contain errors and it is not known which copy of the chromosome a particular fragment belongs to . the problem of haplotyping has been studied extensively . the general minimum error correction ( mec ) problem was proved to be np - hard . it was also proved to be np - hard even if the snp matrix is gapless using a reduction from the max - cut problem . hapcut   and refhap   are two of the most accurate algorithms in this regard . in this paper , we give a heuristic algorithm for individual haplotyping based on minimum error correction . the complexity of each iteration is o(mk ) for an snp matrix of dimension ( m , k ) . the algorithm is inspired from the famous fiduccia and mattheyses ( fm ) algorithm for bipartitioning a hypergraph minimizing the cut size . extensive simulations indicate that hmec outperforms the genetic algorithms of wang et al .   in terms of both reconstruction rate and running time , and it has better ( in most cases ) or comparable accuracy and significantly smaller running time than that of hapcut , which is the most accurate heuristic algorithm available . we also compared hmec with some other algorithms such as speedhap , fasthare , mlf , 2 distance mec , and shr-3   using the hapmap - based instance generator and comparison framework [ 14 , 15 ] . the rest of the paper is organized as follows . in section 2 , we present some definitions and preliminary ideas . in section 3 , we report on an extensive performance study evaluating hmec with other available techniques in section 5 . let f be the set of m fragments produced from two copies of the chromosome . each fragment contains information of nonzero number of snps in s. because the snps are bi - allelic , let the two possible alleles for each snp site be 0 and 1 , where they can be any two elements of the set { a , t , g , c}. since all the nucleotides are the same at the sites other than snp sites , we can remove these extraneous sites from all the fragments and consider the fragments as the sequences of the snp sites only . thus each fragment f  f is a string of symbols { 0,1 ,  } of length k where all the fragments can be arranged in an m  k matrix m = { mij } , i = 1 ,  , m , j = 1,  k , where row i is a fragment from f and column j is a snp from s. this matrix is called the snp matrix as follows 
 ( 1)|110100011101011101001001110100010101101010110101101011| the consecutive sequence of  s that lies between two nonhole symbols is called a gap . a gapless snp matrix is the one that has no gap in any of the fragments . in ( 1 ) , the first , second , and third rows have no gaps while each of the fourth and sixth rows has one gap . , mm can be viewed as an ordered set of m fragments where a fragment mi = mi1 , mi2 ,  , mik is an ordered set of k alleles . a fragment mi is called to cover the jth snp if mij  { 0,1 } and called to skip the jth snp if mij = . let ms and mt be two fragments . the distance between two fragments , d(ms , mt ) , is defined as the number of snps that are covered by both of the fragments and have different alleles . hence , 
 ( 2)d(ms , mt)=j=1kd(msj , mtj ) , 
 where d(x , y ) is defined as 
 ( 3)d(x , y)={1,if  x  and  y  and  xy;0,otherwise . in ( 1 ) , the distance between the second and the third fragment is two , as they differ in the seventh and ninth snp sites ( columns ) . two fragments ms and mt are said to be conflicting if d(ms , mt ) > 0 . let p(c1 , c2 ) be a partition of m , where c1 and c2 are two sets of fragments taken from m so that c1c2 = m and c1c2 =  . in figure 1(b ) , an arbitrary partition corresponding to the snp matrix of figure 1(a ) is shown . a snp matrix m is an error - free matrix if and only if there exists a partition p(c1 , c2 ) of m such that for any two fragments x , y  ci , i  { 1,2 } , x and y are non - conflicting , that is , d(x , y ) = 0 . the partition in the figure 1(b ) is not error free since d(m1 , m2 ) > 0 in c1 and d(m5 , m6 ) > 0 in c2 . for an error - free snp matrix , a haplotype hi , i  { 1,2 } is constructed from its corresponding fragment class ci using the following formula : 
 ( 4)hij={1,if  at  least  one  fragment  in  ci  has  a  1  in  jth  snp;0,if  at  least  one  fragment  in  ci  has  a  0  in  jth  snp;,if  all  the  fragments  in  ci  skips  jth  snp ; 
 where ci is called the defining class of haplotype hi , and hij , where  i  { 1,2 } and j = 1,  k , denotes the jth element of the haplotype hi . if a matrix m is not error - free , there will be no error - free partition p. for such a matrix m , there will be at least one conflicting pair of fragments in each of the classes for all possible partitions . therefore it is impossible to construct a haplotype that is non - conflicting with all the fragments in its defining class of fragments . if we are given a partition p(c1 , c2 ) and two haplotypes h1 and h2 constructed from p then the number of errors e(p ) that needs to be corrected can be calculated by the following formula : 
 ( 5)e(p)=i=12fcid(f , hi ) . the mec problem asks to find a partition p that minimizes the error function e(p ) over all such partitions of an snp matrix m. in this section , we give our heuristic algorithm based on minimum error correction which we call hmec . construction of a haplotype from an erroneous class c requires correction of snp values , that is , alleles , in the fragments . therefore , for each snp site , the haplotype should take the allele that is present in the majority of the fragments . let nj(c ) be the number of fragments of a collection c that have 0 in the jth snp . therefore , to minimize the number of errors e(p ) for a specific partition p , the haplotype should be constructed according to the following methodology : 
 ( 6)hij={1,if  nj1(ci)>nj0(ci);0,if  nj0(ci)nj1(ci )  and  nj0(ci)0;,if  nj1(ci)=nj0(ci)=0 ; 
 where i  { 1,2 } and j = 1,2 ,  , k. in figure 1(c ) , two haplotypes h1 and h2 , associated with the partition p in figure 1(b ) , are constructed in this way . this algorithm starts with a current partition pc = p(m ,  ) and iteratively searches a better partition . in each iteration , the algorithm performs a sequence of transfer of fragments from their present collection to the other one so that the partition becomes less erroneous . the transfer of a fragment from one collection to the other can increase or decrease the error function e(p ) . let the partition before transferring a fragment f be pp and the partition resulted is pn . we define the gain of the transfer as gain(f ) = e(pp )  e(pn ) . , m } be an ordering of all the fragments in a partition p in such a way that fragment fi will precede fragment fj if all the fragments before fi in f have already been transferred to form an intermediate partition pi and gain(fi )  gain(fj ) over pi . we also define the cumulative gain of a fragment ordering f up to the nth fragment as cgain(f , n ) = i=1gain(fi ) . here gain(fi ) = e(pi )  e(pi+1 ) . the maximum cumulative gain , mcgain(f ) , the algorithm finds the current ordering fc of pc and transfers only those fragments of fc that can achieve the mcgain(fc ) and the fragment that is the last to be transferred is referred as fmax . thus the algorithm moves from one partition to another reducing the error function by an amount of mcgain(fc ) . the algorithm continues as long as mcgain(fc ) > 0 and stops whenever mcgain(fc )  0 . first , to find fc in each iteration , the algorithm repeatedly transfers the fragment that is not transferred previously in this iteration and has maximum gain over all such fragments . to accomplish this , we use a locking mechanism . at the beginning of each iteration , all the fragments are set free . the free fragment with maximum gain is found out and tentatively transferred to the other collection . the algorithm then finds the next free fragment with maximum gain in p1 and transfer and lock that fragment to create the p2 . thus , free fragments are transferred until all the fragments are locked and the order of the transfer ( fc ) is stored in the log table along with the cumulative gains ( cgain ) . mcgain is the maximum cgain and fmax is the fragment corresponding to mcgain in the log table . after finishing all such tentative transfers , current  partition of the next iteration , the algorithm checks the log to find the mcgain(fc ) and fmax , and rollback the transfer of all the fragments that were transferred after fmax . while tentatively transferring a free fragment , the algorithm needs to find the fragment with maximum gain among the free fragments ( which are not yet transferred ) . = e(pp )  e(pn ) for a fragment , we need to calculate two error values of two different partitions : the present intermediate partition and the next partition which will be resulted if f is transferred . each of these error function requires calculation of two new haplotypes from their corresponding collections ( see figure 2 ) . although e(pp ) and the haplotypes of pp can be found from the previous transfer , calculation of e(pn ) requires construction of haplotypes of pn . since , the difference between pp and pn is only one transfer , we can introduce differential calculation of haplotypes hi , i  { 1,2 } of next partition from the haplotypes of hi , i  { 1,2 } of present partition . for this purpose , the algorithm stores nj(ci ) and nj(ci ) values of the present partition . after a transfer these values will either remain same or be incremented or decremented by 1 . that is why it is now possible to construct hi , i  { 1,2 } in o(k ) time . to compute e(pn ) from the haplotypes requires o(mk ) time . thus running time to compute the e(pn ) as well as to compute gain(f ) is o(mk + k ) . for each intermediate partition pi , i = 1 , , n , we need to compute gain measures for m  i unlocked fragments to find the maximum one . the transfer of this fragments requires updating of nj(ci ) and nj(ci ) , i  { 1,2 } and j = 1,2 , finally , there will be m such transfer in each iteration and maximum m rollbacks . thus each iteration will require o(m(m(mk + k ) + k ) + mk ) ~ o(mk ) running time . we consider that the current partition pc = p1 is the partition given in figure 1(b ) for the snp matrix m of figure 1(a ) . all the intermediate partitions pi , i  { 1 ,  , 7 } are shown sequentially and the gains of each fragment over the intermediate partitions are shown on the right of each partition . for example , the fragment with the maximum gain in p2 is fragment 6 which has gain two . after each transfer , , the ordering fc of the fragments is 2,6 , 5,1 , 4,3 which is also the order of locking of the fragments . all the tentative transfers after fmax have to be rolled back so that the p3 becomes the next pc . we now give an approximate gain measure to make our algorithm faster . for large snp matrix , we can use an approximation in the calculation of the gain(f ) by using only the fragment f and not using the m  1 other fragments . hi is the haplotype of f 's present collection ci of partition pp , and hj is the haplotype of f 's next collection cj of partition pn . this function ignores the effect of fragments other than f on gain(f ) , but reduces the run time of gain calculation to o(k ) . in this section , we demonstrate the performance of our algorithm using both real biological and simulated datasets . we compared our algorithm with gmec   and hapcut . we performed the simulation using the data from angiotensin - converting enzyme ( ace )   and public daly set   to compare with gmec , and used the huref data   to compare with hapcut . we also compared hmec with some other algorithms ( speedhap , fasthare , mlf , 2 distance mec , shr-3 ) using rehap website interface , which is an hapmap - based instance generator and comparison framework .   in this section , we compare the performance of our algorithm with the genetic algorithm , which we call gmec , described in . we first sample the original haplotype pair into many fragments with different coverage and error rates . coverage rate indicates the percentage of the total columns of the snp matrix that have been sampled out . we varied the error rate while number of fragments and coverage rate were kept constant . also , coverage was varied while number of fragments and error rate were kept constant . notice that the way we introduced error and controlled the coverage rate is not necessarily same as that of . therefore , the reconstruction rates of the branch and bound algorithm described in , which is an exact algorithm , should not be compared with those of hmec . angiotensin - converting enzyme catalyses the conversion of angiotensin i to the physiologically active peptide angiotensin ii , which controls fluid - electrolyte balance and systematic blood pressure . because it has a key function in the renin - angiotensin system , many association studies have been performed with dcp1 ( encode angiotensin - converting anzyme ) . completed the genomic sequencing of the dcp1 gene from 11 individuals and reported 78 snp sites in 22 chromosomes . we generate 50 fragments from each of these haplotype pairs with varying coverage and error rate . we perform the simulation for three different coverage rates ( 25% , 50% , and 75% ) . for each of these coverage rates , we perform our simulation for different error rates such as 5% , 10% , 15% , 20% , 25% , 30% , and 50% . in every case figure 5 illustrates the comparison that bears the clear testimony to the superiority of our algorithm . for most instances , the reconstruction rate achieved by our algorithm is 100% or greater than 98% . only for a few cases with very high error rate and low coverage value we also perform the experiment for different coverage rates while keeping the error rate constant . also , our algorithm clearly outperforms the genetic algorithm except for very low coverage value ( which is unrealistic ) .   in this section , daly et al . reported a high - resolution analysis of a haplotype structure across 500  kb on chromosome 5q31 using 103 snps in a european derived population which consists of 129 trios [ 18 , 20 ] . we performed the experiment exactly in the same way that we did for angiotensin - converting enzyme . experimental results suggest that hmec is much better than the genetic algorithm . again , for most cases , the reconstruction rate achieved by our algorithm is 100% or greater than 98% . for every instance one of the very important advantages of our algorithm is that it takes very short time to reconstruct the haplotypes . , we perform the simulation by varying the length ( length denotes the number of snp sites in the haplotype pair ) of the haplotypes while fixed the value of the coverage rate and the error rate at 50% . since haplotypes with such varying lengths are not available , we rely on the simulated data . clearly , hmec is much faster than gmec . for example , while hmec can reconstruct a haplotype with 936 sites in a fraction of a second , gmec takes 72 seconds . it computes max - cut on the graph to find the position to flip and iterates until no improvement in mec score is achieved . experimental results suggest that although hapcut is reliable , its running time is too large to be a realistic choice for whole genome haplotyping . we used the filtered huref data from levy et al .   to evaluate the performance . we generated several test data sets varying the coverage and error rate , and tested the performance of hmec and hapcut on these data sets . experimental results indicate that the reconstruction rates of both hapcut and hmec are reasonably good . notice that hapcut is only better than hmec for very low ( and thus unrealistic ) coverage values ( 5% and 10% coverage ) . furthermore , hapcut is significantly slower than hmec . for an instance , with 35% coverage and 40% error rate , hapcut takes 215.5 seconds where hmec takes only a fraction of a second ( see table 5 ) . therefore , although generally hapcut provides reliable reconstruction rate , on large dataset , it is an unrealistic choice due to its time consuming operations . on the other hand we created the allele matrix from the rehap error matrix and fed the matrix to both hmec and hapcut algorithms . hmec consistently produces higher reconstruction rate than hapcut ( see table 6 ) . also , the running time of hmec is clearly much better than hapcut . we also compared hmec with some other well - known algorithms such as speedhap , fasthare , mlf , 2 distance mec , and shr-3   using hapmap - based instance generator and comparison framework [ 14 , 15 ] . the results are shown in table 7 and table 8 that suggest that no single method clearly outperforms the others in all cases . however , reconstruction rates achieved by hmec are the highest or very close to the highest . this bears a clear testimony to its suitability as a practical tool for individual haplotyping . in this paper , we present a heuristic algorithm ( hmec ) based on minimum error correction that computes highly accurate haplotypes significantly faster than the known algorithms for haplotyping . the algorithm is inspired from the famous fiduccia and mattheyses ( fm ) algorithm for bipartitioning a hyper graph minimizing the cut size . we report on an extensive performance study evaluating our approach with other available techniques using both real and simulated datasets . comprehensive performance study shows that our algorithm outperforms ( in most cases ) or matches the accuracy of other well - known methods , but runs in a fraction of the time needed for other techniques . high accuracy and very fast running time make our technique suitable for genome - wide scale data . for example , small groups of fragments that are declared to be in the same haplotype can be identified . probabilistic methods like expectation maximization ( em ) also deserve some consideration over such optimization problems . in the near future
haplotype is a pattern of single nucleotide polymorphisms ( snps ) on a single chromosome . 
 constructing a pair of haplotypes from aligned and overlapping but intermixed and erroneous fragments of the chromosomal sequences is a nontrivial problem . 
 minimum error correction approach aims to minimize the number of errors to be corrected so that the pair of haplotypes can be constructed through consensus of the fragments . 
 we give a heuristic algorithm ( hmec ) that searches through alternative solutions using a gain measure and stops whenever no better solution can be achieved . 
 time complexity of each iteration is o(m3k ) for an m  k snp matrix where m and k are the number of fragments ( number of rows ) and number of snp sites ( number of columns ) , respectively , in an snp matrix . 
 alternative gain measure is also given to reduce running time . 
 we have compared our algorithm with other methods in terms of accuracy and running time on both simulated and real data , and our extensive experimental results indicate the superiority of our algorithm over others .