BACKGROUND
linkage analysis is the first step in the search for a disease gene.
the aim is to find the rough location  of the gene in the chromosome.
linkage studies have facilitated the identification of several hundred human genes that can harbor mutations leading to a disease phenotype  <cit> .
the principle of linkage analysis is simple.
all our chromosomes come in pairs, one inherited from the mother and the other from the father.
each pair of chromosomes contains the same genes in the same order, but the sequences are not identical.
thus it is possible to find out whether a particular sequence comes from the mother or father.
these sequence variants are called maternal and paternal alleles.
the key problem for linkage analysis is to infer the pairs of alleles and identify regions whose allele is shared by all or most of the diseased individuals but by none or few of the normal individuals.
linkage analysis has been extensively studied in recent years.
almost all the existing methods are for families with clearly given pedigrees.
the pedigree information helps a lot in the design of computational methods.
early approaches to linkage analysis were based on sparse microsatellite markers.
with the new development of microarray techniques, high-density snp genotype data can be used for large-scale and cost-effective linkage analysis  <cit> .
with high-density snp genotype data, there exists a sufficient number of informative markers between every pair of recombination points, and the allele-sharing status among the family members can be unambiguously determined.
lots of new computer programs have been developed for dealing with high-density snp genotype data.
there are two categories of existing approaches to linkage analysis, the probabilistic approaches and the deterministic approaches.
in probabilistic approaches, recombinant rates are estimated in a way to maximize the likelihood of the observed data  <cit> .
software tools based on this kind of approach include genehunter  <cit> , linkage  <cit> , allegro  <cit> , merlin  <cit> , etc.
according to  <cit> , these tools have different performances and efficiencies.
some of them  do not work well when the number of markers is large, while the others  do not work well with large number of family members.
though tremendous improvement has been made to them through subsequent modifications  <cit> , this still remains a problem in practice.
on the other hand, these tools can give very accurate results when the size of the pedigree is small.
some deterministic approaches have been developed recently.
the main idea is to infer the haplotype segments based on the input genotype data so that all or most of the diseased individuals share a segment that is shared by none of the normal individuals  <cit> .
the mathematical model used here is to minimize the total number of recombinants among all the individuals in the pedigree.
lots of algorithms for haplotype inference with a pedigree have been developed.
qian and beckmann  <cit>  and tapader et al.
<cit>  proposed a method to minimize the number of recombinants with a given pedigree.
zhang et al.
<cit>  developed a program for general pedigrees assuming that there is no recombinant on the segment.
doi et al.
<cit>  designed two algorithms for haplotype inference with a given pedigree.
one of their algorithms works well when the number of marker loci is a fixed constant, while the other works well when the number of family members is bounded by a small constant.
li and jiang  <cit>  proposed to use an integer linear programming approach for minimum recombinant configuration.
xiao et al.
<cit>  designed a faster algorithm for the case where there is no recombinant.
the algorithm in  <cit>  uses a set of heuristics for haplotype inference with a given pedigree and can give very accurate results when the number of family members is large enough and for each nuclear family the genotype data for both parents are available.
an extended software package  was developed in  <cit>  and it focuses on handling the case where the genotype data for the whole chromosome of one of the parents in a nuclear family are missing.
it also uses the minimum recombinant model for haplotype inference in pedigrees.
throughout this paper, we study the dominant inheritance situation, where sharing of one mutation allele can cause a disease phenotype.
we deal with a very important case, where the sampled individuals are closely related, but the pedigree is not given.
this situation happens very often in lots of villages in china when the individuals share a common ancestor  <dig> or more generations ago.
handling this case will be very helpful to identify some local genetic diseases in china.
the situation also happens when studying wild animals, where the pedigree can not be identified.
to our knowledge, no algorithm can give good results for this case.
to solve this problem, we first developed some heuristic algorithms for haplotype inference without any given pedigree.
we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed in  <cit> .
our heuristic algorithm uses clark’s inference rule  <cit>  to infer haplotype segments.
experiments show that our program performs well.
the recall value is from 90% to 99% in various cases.
this implies that the program can report more than 90% of the true mutation regions.
the value of precision varies from 29% to 90%.
when the precision is 29%, the size of the reported regions is three times that of the true mutation region.
this is still very useful for narrowing down the range of the disease gene location.
our program can complete the computation for all the tested cases, where there are about  <dig>  snps on a chromosome, within  <dig> seconds.
implementation
our software is implemented in java.
it takes the genotype data on a chromosome as well as the disease status for a set of input individuals without any pedigree information as input, and outputs the predicted mutation regions.
the software is platform independent.
in the following, we show the methods we used and how we implemented our algorithm.
the problem
suppose that there is a  pedigree containing many  generations, where we only have the genotype data on a chromosome for the individuals in the latest generation .
here those individuals with given genotype data are referred to as the input individuals.
for each input individual, we also know if such an individual is diseased or normal.
an example is given in figure  <dig>  where there are five generations in the pedigree and we only have the genotype data for the individuals in the dashed rectangle at the bottom of the pedigree.
in such a figure, a square represents a male, while a circle represents a female.
moreover, a filled square  represents a diseased male , while an unfilled square  represents a normal male .
furthermore, if two squares  enclose the same number in the figure, then they correspond to the same male  and their sides  are dashed.
for a genotype segment g of length l, the value at each position of g can be  <dig> , or  <dig>  a position of g with  <dig> indicates that both haplotypes have  <dig> at this position, while a position of g with  <dig> indicates that both haplotypes have  <dig> at this position.
if the value at a position is  <dig>  then one of the haplotypes is  <dig> while the other is  <dig> at this position.
a pair of haplotype segments
is a haplotype pair for a genotype segment g if they satisfy the following conditions:
c <dig>  if the value of g is  <dig>  at a position, the values of h and h′
at this position are both  <dig> .
c <dig>  if the value of g is  <dig> at a position, one of h and h′
is  <dig> and the other is  <dig> at this position.
we also say that the pair of haplotype segments
can explaing.
throughout this paper, we study the dominant inheritance situation, where sharing of one mutation allele can cause a disease phenotype.
the general problem is as follows: we are given two sets of genotypes on the whole chromosome d = {g <dig> g <dig> …,gk} and n = {gk +  <dig> gk +  <dig> …,gn}, where the k genotypes in d are from diseased individuals and the n − k
genotypes in n are from normal individuals.
the n individuals in d and n are closely related .
the objective is to detect the mutation regions on the chromosome, where all the diseased individuals share a common haplotype segment on the mutation region and none of the normal individuals has such a common haplotype segment on the mutation region.
note that each individual has two haplotype segments on each region.
if we know the haplotypes of each input individual over the chromosome, the shared mutation regions can be computed by finding the haplotype segments which are shared by all the diseased individuals but by none of the normal ones.
the true mutation region is a shared mutation region containing the disease gene.
therefore, to solve the problem, the key issue is to infer the haplotype segments.
the task of inferring the haplotypes of each individual over the whole chromosome is extremely hard.
for our purpose, we divide the whole chromosome into a set of disjointed length l segments, where l is a parameter to be determined later.
for each length l segment, we try to infer the two haplotype segments of each individual based on the following mathematical model.
mutation region haplotype inference problem
given two sets dr={g1r,g2r,…,gkr} and nr={gk+1r,gk+2r,…,gnr} of genotype segments on a length l region r, where the first k genotype segments in dr
are from diseased individuals and the n − k
genotype segments in nr
are from normal individuals, we want to compute a center haplotype segment hr
and a pair of haplotype segments hi,1r and hi,2r for each gir in dr∪nr, such that the following conditions hold:
hi,1r=hr for any gir∈dr.
ht,rr≠hr for gtr∈nr and r =  <dig> ;
the total number xr
of distinct haplotype segments on r is minimized.
without loss of generality, we assume that all the genotype segments on r in dr are distinct.
similarly, all the genotype segments on r in nr
are also distinct.
however, a genotype segment on r from a diseased individual and a genotype segment on r from a normal individual may be identical.
in this case, such a genotype segment on r should be in both dr and nr.
condition  makes sure that all the diseased individuals have a haplotype segment which is identical to the center haplotype segment hr.
condition  ensures that all the normal individuals do not have hr.
condition  uses the parsimony principle, i.e., we want the total number of distinct haplotype segments to be minimized.
due to condition , this mathematical model can be viewed as an extension of the parsimony model first proposed by gusfield in  <cit>  for haplotype inference.
the parsimony model has been extensively studied in  <cit> .
mrhip can be viewed as a simplified version of the maximum resolution  problem which is proved to be np-hard in  <cit> .
it should be emphasized that for some input of mrhip on a region r, the solution of mrhip may not exist.
even if the solutions exist, the values of xr
may vary for different inputs.
if r is the mutation region, the solution for mrhip on r always exists and the value of xr should be small.
our approach contains three steps.
first, we decompose the whole chromosome into a set of disjointed length l  segments and try to give a solution for mrhip on each length l segment.
we then have an algorithm to merge length l segments based on the computational results to form longer segments and try to get solutions for mrhip on those longer segments.
after that, we have a method to further extend the longer segments to the left and right.
finally, our algorithm reports all the detected mutation regions.
the algorithm for mrhip
for the mutation region haplotype inference problem, we designed an algorithm to solve it.
given an instance of mrhip, there may or may not exist a solution.
if a solution does not exist, there are two cases:
<dig> there does not exist a center haplotype hr
which is shared by all the genotype segments in dr. this case is referred to as type i. type i cases occur when one element in dr
has genotype value  <dig> and the other element in dr
has genotype value  <dig>
<dig> we can find a center haplotype hr, but some genotype segments in nr
must be explained by a pair of haplotype segments and one of the haplotype segments is identical to hr.
this case is referred to as type ii.
an example of a type ii case is the following: dr = {g1 =  <dig> g2 = 121}, and nr = {g3 =  <dig> g4 = 102}.
based on g <dig> and g <dig>  the shared center haplotype hr must be  <dig>  however, g3
indicates that normal individuals also have a haplotype segment  <dig> on r which is identical to the shared center hr.
thus, condition  in mrhip does not hold.
in our algorithm, we first compute the center haplotype hr based on the diseased genotype segments in dr. we look at the positions in r one by one.
based on c <dig> and c <dig>  if one of the diseased individuals has genotype value  <dig>  then the haplotype value of hr
at this position should be 0; if one of the diseased individuals has genotype value  <dig> at a position, then the haplotype value of hr
at this position should be  <dig>  if there exists a position p at which one diseased individual has genotype value  <dig> and the other diseased individual has genotype value  <dig>  then a conflict occurs and position p is called a conflicting position.
once a conflict occurs, we simply conclude that there is no mrhip solution on this segment r. we say the type i false occurs in this case.
if all the diseased individuals have genotype value  <dig> at a position p in r, then the haplotype value of hr
cannot be determined at this step.
we call such a position the wild card position and put a ∗ at the wild card position to indicate that the haplotype value of hr
will be determined later.
the detailed procedure  for computing the center haplotype segment hr is given as follows:
for each position p in rdo
<dig> if all the diseased individuals have genotype value  <dig> or  <dig>  then set the haplotype value of hr
at p to  <dig>
<dig> if all the diseased individuals have genotype value  <dig> or  <dig>  then set the haplotype value of hr
at p to  <dig>
<dig> if some diseased individuals have genotype value  <dig> and some other diseased individuals have genotype value  <dig>  then return type i false.
<dig> if all the diseased individuals have genotype value  <dig>  then set the haplotype value of hr
to ∗
.
without loss of generality, for each diseased individual gir∈dr, we set hi,1r=hr for i ≤ k. then we can set hi,2r in such a way that  is a haplotype pair for gir with i ≤ k. note that, the values of  at wild card positions in r are still not yet determined.
here we refer to hi,2r as a partially inferred haplotype segment on r. let intq={h <dig> r,h <dig> r,…,hk,2r}, where if two haplotype segments hi,2r and hi′,2r are identical, then we just keep one of them.
note that if any haplotype segment hi,2r∈intq is undetermined at a position p, then all the haplotype segments in intq are undetermined at p.
after partially determining hr
and  for every gir in dr, we use a heuristic method to infer the haplotype segments for gir∈nr.
since we want to minimize the total number of resulting distinct haplotype segments on r, our strategy is to let the inferred haplotype segments for gir∈nr share as many haplotype segments as possible.
this is actually clark’s inference rule  <cit> .
let q be a queue that contains a set of  inferred haplotype segments on r. initially, q = intq.
a partially inferred haplotype segment h in q can solvegjr if the following conditions hold:
<dig> if h is  <dig> at a position p then gjr is  <dig> or  <dig> at position p.
<dig> if h is  <dig> at a position p then gjr is  <dig> or  <dig> at position p.
we can use h to solvegjr by constructing two haplotype segments  as follows:
using h to solve gjr:
<dig> if h is  <dig> at position p in r then we set hj,1r= <dig> at p and hj,2r at p is set according to rules c <dig> and c <dig>
<dig> if h is  <dig> at position p in r then we set hj,1r= <dig> at p and hj,2r at p is set according to rules c <dig> and c <dig>
<dig> if h is undetermined at position p in r, and gjr is  <dig>  at p, then set hj,1r=h=hj,2r= <dig>  at p. here we also have to determine the value of h at position p accordingly.
after the undetermined value of h at p is determined, if h is obtained from a gir∈dr, then we also have to determine the values of hr
and other hi,2r for gir∈dr at position p according to the haplotype value of h at p and rules c <dig> and c <dig>
<dig> if h is undetermined at position p in r, and gjr is  <dig> at p, then hj,1r and hj,2r remain undetermined at p.
a genotype segment gir∈nr is solved if the pair of haplotype segments  for gir are  determined.
in our algorithm, we use p to store the set of genotypes in nr that have not been solved.
initially, p = nr.
we then use the haplotype segments in q one by one and try to solve each of the genotypes in p. after trying to use a h ∈ q to solve all gjr’s in p, we delete h from q. two haplotype segments on r in q are compatible if there does not exist any position p in r such that one segment has value  <dig> and the other segment has value  <dig> at p. for two compatible haplotype segments h <dig> and h <dig> on r, we can merge them to form one haplotype segment h, where the value of h is determined as  <dig> or  <dig> if at least one of h <dig> and h <dig> is  <dig> or  <dig> and the value of h remains undetermined if both h1
and h <dig> are undetermined.
again, for any position p, if the value of h1
or h <dig> is not identical to that of h, then the value of h1
or h <dig> is changed from undetermined to  <dig> or  <dig>  thus, we have to update some of the previously inferred haplotype segments accordingly.
when we use a h ∈ q to solve a gjr in p, we can obtain another new haplotype segment h′
on r. if h′
is compatible with a haplotype segment in q, we then merge them.
note that, h′
might be compatible with more than one haplotype segment in q. in this case, we arbitrarily choose a compatible haplotype segment in q and merge the two haplotype segments.
if h′ is not compatible with any haplotype segment in q, we add h′
into q.
we give an example to illustrate the above process.
example 1
dr = {g1 =  <dig> g2 =  <dig> g3 = 10221}
and nr = {g4=10121}.
after procedure p <dig>  hr = 10∗11
and consequently intq = {h <dig>  = 10∗ <dig> h <dig>  = 11∗ <dig> h <dig>  = 10∗01}.
after that , we can use h <dig>  = 10∗11
in intq to solve g4 = 10121
in nr.
based on h <dig>  = 10∗ <dig>  g4 = 10121
can be solved as h <dig>  = 10111
and h <dig>  =  <dig>  moreover, since we want h <dig>
and h <dig>
to be identical , h <dig>
is updated as h <dig>  =  <dig>  correspondingly, we update hr = h <dig>  =  <dig>  h <dig>  =  <dig>  h <dig>  =  <dig> and h <dig>  =  <dig>  after that, the set of distinct haplotype segments we have obtained so far is {hr =  <dig> h <dig>  = h <dig>  =  <dig> h <dig>  =  <dig> h <dig>  =  <dig> h <dig>  = 10101}.
note that h <dig>
and h <dig>
are compatible , and the set of distinct haplotype segments is {hr =  <dig> h <dig>  = h <dig>  =  <dig> h <dig>  =  <dig> h <dig>  = h <dig>  = 10101}.
after trying to use all the h’s in q to solve all gjr’s in p, q will become empty.
when q is empty and p still contains at least two genotype segments, we consider all pairs of genotype segments gjr and gj′r in p and use the following method to infer the haplotype segments.
inferring the haplotype segments from a pair gjr and gj′r in p:
a position p in r is a conflicting position for gjr and gj′r if one of gjr and gj′r has the genotype value  <dig> and the other has genotype value  <dig> at p. the pair of gjr and gj′r can share a common haplotype segment on r if there is no conflicting position in r for gjr and gj′r.
the shared haplotype segment can be computed as follows:  if one of the genotype values at position p is  <dig>  then the haplotype value is  <dig> at p;  if one of the genotype values is  <dig> at p, then the haplotype value is  <dig> at p; if both genotype values are  <dig> at p, then the haplotype value at p is undetermined.
once the shared haplotype segment for gjr and gj′r are computed, we can determine the other haplotype segments for gjr and gj′r based on c <dig> and c <dig>
after inferring the haplotype segments from a pair gjr and gj′r that can share a common haplotype segment, we delete gjr and gj′r from p, merge compatible inferred haplotype segments, and insert the newly obtained haplotype segments into q. once q is not empty, we can use haplotype segments in q to solve the genotype segments in p again.
the process is repeated until p is empty.
the detailed algorithm is given as algorithm  <dig>
algorithm 1
mutation region haplotype inference
input: two sets of genotype segments dr={g1r,g2r,…,gkr} and nr={gk+1r,gk+2r,…,gnr} on r.
output:true if there is a solution.
type i false if two diseased haplotype segments are conflict at a position in r; type ii false otherwise.
1: compute the center haplotype hr
as in procedure p <dig>
2:
ifhr
does not exist then
3:
return type i false
4:
else
5:
set hi,1r=hr for gir∈dr;
6: compute hi,2r according to hi,1r and c <dig> and c <dig> for each gir∈dr; set q={h <dig> r,h <dig> r,…,hk,2r}  and p = nr.
7:
end if
8:
whileq ≠ ∅
and p ≠ ∅do
9: delete a haplotype segment hi,2r from q;
10:
ifhi,2r can solve gjrthen
11: use hi,2r to solve gjr.
add the newly obtained haplotype segments  into q and delete gjr from p.
12:
end if
13:
end while
14:
if there are at least  <dig> genotype segments in pthen
15:
if there exists a pair of genotype segments gjr and gj′r that can share a haplotype segment on rthen
16:
infer the haplotype segments of gjr and gj′r and insert them  into q, goto line  <dig>
17:
end if
18:
if any inferred haplotype segments for some gjr∈nr on r is identical to hrthen
19:
return type ii false.
20:
end if
21:
fix hs,1r and hs,2r for each gsr∈p so that hs,1r≠hr and hs,2r≠hr.
22:
if line line  <dig> fails then
23:
return type ii false
24:
else
25:
return true.
26:
end if
27:
end if
the following is an example to illustrate the case when q becomes empty.
example 2
dr = {g1 =  <dig> g2 = 11210}
and nr = {g3 =  <dig> g4 = 11121}.
after procedure p <dig>  hr = 11110
and intq = { <dig> }.
after trying all the h’s in intq to solve gi’s in nr, q becomes empty and p = nr  ={ <dig> }.
in this case, we look at both  <dig> and  <dig> in p and infer a shared haplotype h <dig>  = h <dig>  = 11111
and the other two haplotype segments h <dig>  = 01111
and h <dig>  =  <dig>
the algorithm for the whole chromosome
for an input segment r on a chromosome, if algorithm  <dig> returns true, then r is a valid segment.
in order to get the mutation regions, we decompose the whole chromosome into a set of disjointed length l segments.
for each segment, we run algorithm  <dig> to test if the segment is valid.
after finding all the valid segments, we repeatedly merge two valid segments into a long valid segment if the two segments are within 2l snps and algorithm  <dig> returns type ii false on all the segments in the gap.after the above merging process, we obtain several long valid segments.
for each such long valid segment [sb,se), we run algorithm  <dig> on the three segments [sb− <dig> l,se +  <dig> l), [sb− <dig> l,se +  <dig> l) and [sb,se) and select the longest one ) which returns true.
since we impose that algorithm  <dig> returns type ii false for the segments in gaps in the merging process, we can always ensure that algorithm  <dig> returns true for [sb,se).
extending [b,e) to the left and right:
after we obtain r = [b,e) as discussed above, we try to extend the segment [b,e) to the left and right.
on the segment r=[b,e), we have inferred hi,1r and hi,2r for each gir∈dr∪nr.
these hi,1r and hi,2r form a collection of disjointed sets h <dig> h <dig> …,hm, where each hk is a set of identical haplotypes in {hi,1r,hi,2r|i= <dig> ,…,n} on r.
we extend the segment [b,e)
to the left and right by looking at each position p. we first try p = b − q
for q =  <dig> ,…
and then p = e −  <dig> + q
for q =  <dig> ,…
.
for each hk,
<dig> if there exist some hi,j
in hk
such that gi
are  <dig>  and others are  <dig>  then every hi,j
in hk
should be  <dig> .
if there exists a hi,j
in hk
that has been set to the conflict value  <dig>  before, then we know that position p is a conflicting position and p should not be extended to be part of r and the extension process to the current direction  should stop.
otherwise, we set every hi,j
in hk
to be  <dig> .
if p is not a conflicting position, after setting hi,j
in hk
to be  <dig> , we can determine the value of hi,j′  according to the value of hi,j
at p, the value of gi
at p and rules c1
and c <dig>  again, we should test if such a value of hi,j′ is consistent with the value of hi,j′ determined before .
if conflict exists, then p is a conflicting position and the extension process to the current direction  should stop.
let hi,j′∈hk′.
if there is no conflict, we should also update the value of all h’s in hk′.
this recursive process continues until no further change can be made.
<dig> if all hi,j
in hk
are  <dig>  we set hi,j
as undetermined.
the extension process stops when we find conflicts in both directions.
the extended region obtained from .
after the extension process, our program reports all the mutation regions obtained in the algorithm.
the complete algorithm to find the mutation regions on the whole chromosome is shown as algorithm  <dig>
algorithm 2
the algorithm for the whole chromosome
input: two sets of genotype on the whole chromosome d = {g <dig> g <dig> …,gk}
and n = {gk +  <dig> gk +  <dig> …,gn}.
output: the detected mutation regions
1:decompose the whole chromosome into segments of length l =  <dig>
2:for each segment do
3:
use algorithm  <dig> to test if the segment is valid.
4:end for
5:merge the valid segments  to form longer segments.
6:for each segment [sb,se)
obtained in line  <dig> do
7:
select the longest segment of [sb −  <dig> l,se +  <dig> l), [sb −  <dig> l,se +  <dig> l)
and [sb,se)
which will return true by calling algorithm  <dig>  denote it as [b,e).
8:
extend [b,e)
to get a candidate mutation region.
9:end for
10:output all the mutation regions obtained in the for loop of lines 6- <dig>
RESULTS
in this section, we first show some experiments on simulated data.
we then give a real case study to show that our program can also handle real data .
a discussion is given at the end of this section.
experiments on simulated data
in order to evaluate the performance of our method and the feasibility of the mathematical model proposed in this paper, we write a program in c++ to produce simulated data.
the program takes a pedigree  and the haplotype data for the whole chromosome of each founder in the pedigree as the input.
it generates the haplotype data for the remaining individuals in the pedigree using the standard χ <dig> model for recombination with m  equal to  <dig>  and according to the male/female averaged genetic map for chromosome  <dig> downloaded from hapmap .
also see  <cit> .
the haplotype data of a non-founder in the pedigree are generated to randomly inherit one strand of the four-strand chromatid bundle from each parent of the non-founder.
a mutation point is selected uniformly at random from the snp sites of the chromosome.
each diseased offspring is forced to inherit  the strand with the mutation point and the normal offspring are forced to inherit the strand without the mutation point.
in this way, we can guarantee that there is exactly one true mutation region.
note that the true mutation region must be a shared mutation region.
see implementation for the definition.
moreover, since we know the haplotype data of all the individuals in the simulations, we can easily find the shared mutation regions.
by definition, there may exist more than one shared mutation region.
to generate the simulated data, we randomly chose some of the haplotype data for chromosome  <dig> of  <dig> unrelated japanese in tokyo and han chinese in beijing in the database of hapmap project  as the haplotype data for each founder.
there are about  <dig>  snp sites on this chromosome.
recall that our program takes two sets of individuals d and n and their genotype data as input.
after generating the haplotype data of each individual, we only use some of the individuals  and their associated genotype  data as the input of our program.
to evaluate the performance of our method, we used different pedigrees to evaluate our algorithm.
figures  <dig>   <dig>   <dig> and  <dig> are pedigrees of  <dig> generations with  <dig> to  <dig> diseased individuals in the dashed rectangle at the bottom.
those individuals in the dashed rectangle at the bottom of each pedigree are the input of our program.
the correctly detected mutation regions are the intersection of the regions reported by the computer program and the true mutation region.
here, precision is defined as the number of snps in the correctly detected mutation regions divided by the total number of snps in the regions output by the program.
the value of recall is defined as the number of snps in the correctly detected mutation regions divided by the total number of snps in the true mutation region.
so, if the value of recall is  <dig>  then all the snps in the true mutation region have been reported by the program.
similarly, if precision is  <dig>  then all the reported snps are in the true mutation region.
we performed  <dig> experiments for each pedigree.
since there are about  <dig>  snp sites on the chromosome, we set l =  <dig>  for each region
reported by our algorithm, we define a score as follows: let dh be the number of distinct haplotypes on this region and length =
the length of this region.
then the score of this region is defined as score = ∗length, where n is the total number of input genotypes in d∪n.
this score can balance the length of the mutation region and the number of distinct haplotype segments on the region.
with the longer region and smaller dh, the score becomes higher.
to illustrate the quality of our program, we report the results when our program reports the region with the highest score and the first three regions with the highest scores, respectively.
in fact, our program does not need this score in the computation.
the program simply reports all the mutation regions.
see the genotype data error handling in discussion.
the precision and recall on the experiments are shown in table  <dig>  only the genotype of the individuals in the latest generation of pedigree 1−4
are known in this experiment.
several mutation regions may be detected by our algorithm.
in table  <dig>  the results when our program reports the region with the highest score are shown in the columns under “one region”.
the results when our program reports three regions with the highest scores are shown in the columns under “three regions”.
the precision and recall are calculated based on the true mutation region, the reported region, and the intersection of the reported region and the true mutation region.
the precision’ and recall’ are calculated by replacing the true mutation region with shared mutation regions.
the column “time” indicates the average time of our program by running  <dig> experiments on each pedigree.
from table  <dig>  we can see that the values of recall are from  <dig> % to  <dig> % and the values of precision are from  <dig> % to  <dig> % in the four pedigrees.
when the number of diseased individuals is increased, the values of both recall and precision are improved significantly.
when there are  <dig> or  <dig> diseased individuals, the value of recall is more than 97%.
that is, the program can report most of the snps in the true mutation region.
in practice, one can often get the genotype data for the individuals of the latest two generations.
thus, we study this case by looking at different input individuals based on the pedigrees in figures  <dig>   <dig>   <dig> and  <dig>
now, we study different sets of input individuals in the latest two generations of pedigree  <dig>  these different sets of input individuals in the latest two generations in the pedigree are given in figure  <dig>  a square/circle with a slash indicates that such individual is not included as part of the input though the individual is used in generating the simulated data.
for the rest of test, we performed  <dig> experiments for each case and show the average values.
table  <dig> shows the results for the different sets of input individuals in figure  <dig>  the individuals in the latest and latest two generations are not distinguished in our algorithm.
we just input the genotype for all the individuals without the slash.
again, the setting is similar to that of table  <dig>  from table  <dig>  we can see that the values of recall for different inputs are close to 99% except for 2d-3fam- <dig>  2d-2fam- <dig> and 2d-2fam- <dig>  where the input contains only  <dig> or  <dig> diseased individuals.
comparing table  <dig> and table  <dig>  we can see that more input individuals do help improve the values of precision and recall.
we also performed similar experiments for pedigree 2- <dig> .
the results are similar to that in table  <dig> and are given in additional file  <dig>
we also tested the program using pedigrees containing  <dig> and  <dig> generations and  <dig>   <dig>   <dig>   <dig> diseased individuals, respectively, in the latest generation.
the four pedigrees containing  <dig> generations are shown in additional file 1: figure s <dig>  figure s <dig>  figure s <dig> and figure s <dig> in the additional file.
the four pedigrees containing  <dig> generations are shown in additional file 1: figure s <dig>  figure s <dig>  figure s <dig> and figure s <dig> in the additional file.
again, the input individuals are the individuals in the dashed rectangle at the bottom of the pedigree.
the experiment results for  <dig> generations and  <dig> generations are shown in table  <dig> and table  <dig>  respectively.
the settings of table  <dig> and table  <dig> are similar to that of table  <dig>  table  <dig> and table  <dig> show that the performance of our program for  <dig> and  <dig> generations is similar  to that for  <dig> generations.
similar to the case of  <dig> generations, for pedigrees with  <dig> and  <dig> generations, we also tested various cases when some individuals of the latest two generations are available as input individuals.
the results for  <dig> and  <dig> generations are similar to that of  <dig> generations.
the detailed results are given in the additional file.
a real case study
to illustrate the usefulness of our program, we applied our method to a set of real data originally from the phase ii hapmap database and was studied in  <cit> .
in  <cit> , the authors studied two ceu  families  ceph  <dig> and ceph  <dig> .
they identified a segment  on chromosome  <dig> shared by the four individuals na <dig>  na <dig>  na <dig> and na <dig>  for this set of data, there are totally  <dig>  snp sites on the chromosome after the unknown genotypes are eliminated from the database.
there are totally  <dig>  snp sites between 107m and 110m on chromosome  <dig> starting at the 122348-th snp site and ending at the 128866-th snp site.
we applied our program with the six individuals in the two families ceph  <dig> and ceph  <dig> as input and set the four individuals na <dig>  na <dig>  na <dig> and na <dig> as diseased individuals.
we set l =  <dig>  our algorithm found several segments shared by the four diseased individuals.
the lengths of all the reported segments are approximately  <dig> snp sites except the longest ones.
all these segments are shown in table  <dig>  table  <dig> shows the starting and ending point of the segments, the number of distinct haplotypes on the segment , and the score for each segment.
as there is only one conflicting position  <dig> between segments  <cit>  and  <cit> , we should consider such a conflicting position as a data error.
therefore, segment  <cit>  should be the predicted mutation.
this segment starts at the 124561-th snp site and ends at the 129451-th snp site.
the details are shown in figure  <dig>  we can see that the shared segment found by plink in  <cit>   starts at the 122348-th snp site while the starting position of our reported segment  is  <dig>  for the subsegment  <cit>  that we did not report, we found  <dig> conflicting positions in this subsegment containing  <dig>  snp sites.
however, on the segment  <cit>  containing  <dig>  snp sites reported by our program, there is only one conflicting position.
this is strong evidence that the subsegment  <cit>  is not shared by all the four diseased individuals.
we also looked at the segment  <cit>  with length of  <dig>  snp sites on the right of our reported segment  <cit> , and found  <dig> conflicting positions among the  <dig>  snp sites.
we can see that the subsegments  <cit>  and  <cit>  on the left and right of our reported segment have approximately the same number of conflicting positions.
discussion
haplotype inference methods
as discussed before, if the haplotype data for each input individual are known, the problem of finding the true mutation region is straightforward.
currently, there are several population-based phasing methods that can give accurate haplotype segments  <cit> .
however, these methods can only phase a small number of snps effectively and take an extremely long time to infer the haplotype for the whole chromosome.
lrp in  <cit>  can phase more than  <dig>  snps simultaneously within a reasonable time.
however, it is still very slow for phasing  <dig>  snps of a whole chromosome .
moreover, they cannot directly report the true mutation region for a set of input individuals.
on the other hand, our program can complete the computation in less than  <dig> seconds for about  <dig>  snps with about  <dig> to  <dig> individuals.
related mutation region detection methods
to our knowledge, all the existing software packages  need a clearly given pedigree as part of the input.
if the pedigree is not known, most of the software packages do not work.
our algorithm deals with the case where the input individuals are closely related but the pedigree is not given.
merlin is widely used for linkage analysis, where a pedigree is required as part of the input.
it works well on snp data due to the use of sparse trees.
however, it can only analyze pedigrees of moderate size.
when the family size is big, a large memory space is needed and the computation cannot successfully be completed.
as shown in  <cit> , merlin cannot report the results for some pedigrees, e.g., p <dig> and p <dig> in  <cit> , where there are less than  <dig> input individuals.
however, our program can deal with the cases where the number of input individuals is large.
the first row in table  <dig> of the additional file shows the results for  <dig>  input individuals.
we can see that our method can give very high precision and recall in this case .
therefore, our algorithm can handle some cases which cannot be handled by merlin.
the rule-based algorithm in  <cit>  uses a set of heuristics for haplotype inference with a given pedigree.
it can give very accurate results when the number of family members is large enough and for each nuclear family the genotype data for both parents are available.
however, it does not work well when the genotype data of one of the parents are missing in the nuclear family.
if the data for both parents are missing, it does not work.
liden  <cit>  is an extended software package of the algorithm in  <cit> .
it focuses on handling the case where the genotype data for one of the parents in a nuclear family are missing for the entire chromosome.
but it still does not work when the genotype data for both parents in a nuclear family are missing or the family pedigree is not given in the input.
plink in  <cit>  and beagle in  <cit>  can identify the shared haplotype segment between two individuals based on population-based linkage analysis.
but it cannot automatically identify the mutation region taking a set of individuals as the input, which is expected to give more precise prediction.
the above real case study has illustrated this.
genotype data error handling
the real datasets often contain errors.
handling the genotype data errors is an important issue in practice.
for our program, we have a pre-process step to delete all the snps containing missing data.
that is, if the genotype data for an input individual at an snp site are missing, we delete this snp site from the input.
without this step, we cannot get reasonable results for the real case study.
when the genotype data contain errors, it is hard to detect and correct them.
the errors may affect our program’s results in two ways:  an snp site in the true mutation region may become a conflicting position due to error; and  the number of distinct haplotype segments to explain the genotype data is increased.
when  occurs, our score for the detected mutation regions becomes worse.
when  occurs once, our program reports two detected regions with the conflicting position in between.
see  <cit>  and  <cit>  in table  <dig> of the real case study.
when this kind of error occurs many times, our program reports many regions separated by a few snps in the middle.
when the user looks at the results of our program, it is possible to realize that the few snps between two closely located reported regions are due to errors.
this is similar to other linkage analysis programs such as merlin, where each snp site has a score, and the user decides a region  with high scores as the true mutation region.
our program may work for the situation where the input individuals are from multiple families.
our algorithm tries to find regions shared by all the diseased individuals.
thus, as long as the diseased individuals from multiple families share the same  haplotype segment on the true mutation region, our program should be able to find such region.
even if the haplotype segments from different families on the true mutation region are slightly different, the program should be able to report several smaller regions with a few missing snps in the middle.
again, it is possible for users to figure out the whole true mutation by ignoring the few missing snps in the middle.
CONCLUSIONS
we developed a software package for linkage analysis where the input individuals are closely related, but the pedigree is not known.
we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed by dan gusfield .
our heuristic algorithm simply uses clark’s inference rule to infer haplotype segments.
experiments show that our program can give very high value  of recall in various cases.
this implies that the program can report more than 90% of the true mutation region.
the value of precision varies from 29% to 90%.
when the precision is 29%, the size of the reported regions is three times that of the true mutation region.
this is still very useful for narrowing down the range of the disease gene.
availability and requirements
project name: mrd
project homepage:http://www.cs.cityu.edu.hk/~wenjuacui/software/mutationregion/index.html.
the source code is also available.
operating system: platform independent
programming language: java
other requirements: java  <dig> .
<dig> or higher
license: none
any restrictions to use by non-academics: none
competing interests
the authors declare that they have no competing interests
author’s contributions
lw proposed the topic and ideas for algorithms, wc implemented the program, and both authors devised and developed the method and prepared the manuscript.
both authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary material: this file includes several figures and additional experimental results mentioned in the paper.
it contains the different set of input individuals on pedigree 2- <dig> in the paper, the pedigrees containing  <dig> and  <dig> generations and  <dig> , <dig>  diseased individuals in the latest generation respectively.
the tables show the results on the input of the above figures.
click here for file
acknowledgements
this work is supported by a grant from the research grants council of the hong kong special administrative region, china .
lusheng wang is the corresponding author.

 linkage analysis is the first step in the search for a disease gene. the aim is to find the rough location of the gene in the chromosome. linkage studies have facilitated the identification of several hundred human genes that can harbor mutations leading to a disease phenotype. the principle of linkage analysis is simple. all our chromosomes come in pairs, one inherited from the mother and the other from the father. each pair of chromosomes contains the same genes in the same order, but the sequences are not identical. thus it is possible to find out whether a particular sequence comes from the mother or father. these sequence variants are called maternal and paternal alleles. the key problem for linkage analysis is to infer the pairs of alleles and identify regions whose allele is shared by all or most of the diseased individuals but by none or few of the normal individuals. linkage analysis has been extensively studied in recent years. almost all the existing methods are for families with clearly given pedigrees. the pedigree information helps a lot in the design of computational methods. early approaches to linkage analysis were based on sparse microsatellite markers. with the new development of microarray techniques, high-density snp genotype data can be used for large-scale and cost-effective linkage analysis. with high-density snp genotype data, there exists a sufficient number of informative markers between every pair of recombination points, and the allele-sharing status among the family members can be unambiguously determined. lots of new computer programs have been developed for dealing with high-density snp genotype data. there are two categories of existing approaches to linkage analysis, the probabilistic approaches and the deterministic approaches. in probabilistic approaches, recombinant rates are estimated in a way to maximize the likelihood of the observed data. software tools based on this kind of approach include genehunter, linkage, allegro, merlin, etc. according to, these tools have different performances and efficiencies. some of them do not work well when the number of markers is large, while the others do not work well with large number of family members. though tremendous improvement has been made to them through subsequent modifications, this still remains a problem in practice. on the other hand, these tools can give very accurate results when the size of the pedigree is small. some deterministic approaches have been developed recently. the main idea is to infer the haplotype segments based on the input genotype data so that all or most of the diseased individuals share a segment that is shared by none of the normal individuals. the mathematical model used here is to minimize the total number of recombinants among all the individuals in the pedigree. lots of algorithms for haplotype inference with a pedigree have been developed. qian and beckmann and tapader et al. proposed a method to minimize the number of recombinants with a given pedigree. zhang et al. developed a program for general pedigrees assuming that there is no recombinant on the segment. doi et al. designed two algorithms for haplotype inference with a given pedigree. one of their algorithms works well when the number of marker loci is a fixed constant, while the other works well when the number of family members is bounded by a small constant. li and jiang proposed to use an integer linear programming approach for minimum recombinant configuration. xiao et al. designed a faster algorithm for the case where there is no recombinant. the algorithm in uses a set of heuristics for haplotype inference with a given pedigree and can give very accurate results when the number of family members is large enough and for each nuclear family the genotype data for both parents are available. an extended software package was developed in and it focuses on handling the case where the genotype data for the whole chromosome of one of the parents in a nuclear family are missing. it also uses the minimum recombinant model for haplotype inference in pedigrees. throughout this paper, we study the dominant inheritance situation, where sharing of one mutation allele can cause a disease phenotype. we deal with a very important case, where the sampled individuals are closely related, but the pedigree is not given. this situation happens very often in lots of villages in china when the individuals share a common ancestor or more generations ago. handling this case will be very helpful to identify some local genetic diseases in china. the situation also happens when studying wild animals, where the pedigree can not be identified. to our knowledge, no algorithm can give good results for this case. to solve this problem, we first developed some heuristic algorithms for haplotype inference without any given pedigree. we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed in. our heuristic algorithm uses clark’s inference rule to infer haplotype segments. experiments show that our program performs well. the recall value is from 90% to 99% in various cases. this implies that the program can report more than 90% of the true mutation regions. the value of precision varies from 29% to 90%. when the precision is 29%, the size of the reported regions is three times that of the true mutation region. this is still very useful for narrowing down the range of the disease gene location. our program can complete the computation for all the tested cases, where there are about snps on a chromosome, within seconds. implementation our software is implemented in java. it takes the genotype data on a chromosome as well as the disease status for a set of input individuals without any pedigree information as input, and outputs the predicted mutation regions. the software is platform independent. in the following, we show the methods we used and how we implemented our algorithm. the problem suppose that there is a pedigree containing many generations, where we only have the genotype data on a chromosome for the individuals in the latest generation. here those individuals with given genotype data are referred to as the input individuals. for each input individual, we also know if such an individual is diseased or normal. an example is given in figure where there are five generations in the pedigree and we only have the genotype data for the individuals in the dashed rectangle at the bottom of the pedigree. in such a figure, a square represents a male, while a circle represents a female. moreover, a filled square represents a diseased male, while an unfilled square represents a normal male. furthermore, if two squares enclose the same number in the figure, then they correspond to the same male and their sides are dashed. for a genotype segment g of length l, the value at each position of g can be, or a position of g with indicates that both haplotypes have at this position, while a position of g with indicates that both haplotypes have at this position. if the value at a position is then one of the haplotypes is while the other is at this position. a pair of haplotype segments is a haplotype pair for a genotype segment g if they satisfy the following conditions: c if the value of g is at a position, the values of h and h′ at this position are both. c if the value of g is at a position, one of h and h′ is and the other is at this position. we also say that the pair of haplotype segments can explaing. throughout this paper, we study the dominant inheritance situation, where sharing of one mutation allele can cause a disease phenotype. the general problem is as follows: we are given two sets of genotypes on the whole chromosome d = {g g …,gk} and n = {gk + gk + …,gn}, where the k genotypes in d are from diseased individuals and the n − k genotypes in n are from normal individuals. the n individuals in d and n are closely related. the objective is to detect the mutation regions on the chromosome, where all the diseased individuals share a common haplotype segment on the mutation region and none of the normal individuals has such a common haplotype segment on the mutation region. note that each individual has two haplotype segments on each region. if we know the haplotypes of each input individual over the chromosome, the shared mutation regions can be computed by finding the haplotype segments which are shared by all the diseased individuals but by none of the normal ones. the true mutation region is a shared mutation region containing the disease gene. therefore, to solve the problem, the key issue is to infer the haplotype segments. the task of inferring the haplotypes of each individual over the whole chromosome is extremely hard. for our purpose, we divide the whole chromosome into a set of disjointed length l segments, where l is a parameter to be determined later. for each length l segment, we try to infer the two haplotype segments of each individual based on the following mathematical model. mutation region haplotype inference problem given two sets dr={g1r,g2r,…,gkr} and nr={gk+1r,gk+2r,…,gnr} of genotype segments on a length l region r, where the first k genotype segments in dr are from diseased individuals and the n − k genotype segments in nr are from normal individuals, we want to compute a center haplotype segment hr and a pair of haplotype segments hi,1r and hi,2r for each gir in dr∪nr, such that the following conditions hold: hi,1r=hr for any gir∈dr. ht,rr≠hr for gtr∈nr and r = ; the total number xr of distinct haplotype segments on r is minimized. without loss of generality, we assume that all the genotype segments on r in dr are distinct. similarly, all the genotype segments on r in nr are also distinct. however, a genotype segment on r from a diseased individual and a genotype segment on r from a normal individual may be identical. in this case, such a genotype segment on r should be in both dr and nr. condition makes sure that all the diseased individuals have a haplotype segment which is identical to the center haplotype segment hr. condition ensures that all the normal individuals do not have hr. condition uses the parsimony principle, i.e., we want the total number of distinct haplotype segments to be minimized. due to condition, this mathematical model can be viewed as an extension of the parsimony model first proposed by gusfield in for haplotype inference. the parsimony model has been extensively studied in. mrhip can be viewed as a simplified version of the maximum resolution problem which is proved to be np-hard in. it should be emphasized that for some input of mrhip on a region r, the solution of mrhip may not exist. even if the solutions exist, the values of xr may vary for different inputs. if r is the mutation region, the solution for mrhip on r always exists and the value of xr should be small. our approach contains three steps. first, we decompose the whole chromosome into a set of disjointed length l segments and try to give a solution for mrhip on each length l segment. we then have an algorithm to merge length l segments based on the computational results to form longer segments and try to get solutions for mrhip on those longer segments. after that, we have a method to further extend the longer segments to the left and right. finally, our algorithm reports all the detected mutation regions. the algorithm for mrhip for the mutation region haplotype inference problem, we designed an algorithm to solve it. given an instance of mrhip, there may or may not exist a solution. if a solution does not exist, there are two cases: there does not exist a center haplotype hr which is shared by all the genotype segments in dr. this case is referred to as type i. type i cases occur when one element in dr has genotype value and the other element in dr has genotype value we can find a center haplotype hr, but some genotype segments in nr must be explained by a pair of haplotype segments and one of the haplotype segments is identical to hr. this case is referred to as type ii. an example of a type ii case is the following: dr = {g1 =  g2 = 121}, and nr = {g3 =  g4 = 102}. based on g and g the shared center haplotype hr must be however, g3 indicates that normal individuals also have a haplotype segment on r which is identical to the shared center hr. thus, condition in mrhip does not hold. in our algorithm, we first compute the center haplotype hr based on the diseased genotype segments in dr. we look at the positions in r one by one. based on c and c if one of the diseased individuals has genotype value then the haplotype value of hr at this position should be 0; if one of the diseased individuals has genotype value at a position, then the haplotype value of hr at this position should be if there exists a position p at which one diseased individual has genotype value and the other diseased individual has genotype value then a conflict occurs and position p is called a conflicting position. once a conflict occurs, we simply conclude that there is no mrhip solution on this segment r. we say the type i false occurs in this case. if all the diseased individuals have genotype value at a position p in r, then the haplotype value of hr cannot be determined at this step. we call such a position the wild card position and put a ∗ at the wild card position to indicate that the haplotype value of hr will be determined later. the detailed procedure for computing the center haplotype segment hr is given as follows: for each position p in rdo if all the diseased individuals have genotype value or then set the haplotype value of hr at p to if all the diseased individuals have genotype value or then set the haplotype value of hr at p to if some diseased individuals have genotype value and some other diseased individuals have genotype value then return type i false. if all the diseased individuals have genotype value then set the haplotype value of hr to ∗. without loss of generality, for each diseased individual gir∈dr, we set hi,1r=hr for i ≤ k. then we can set hi,2r in such a way that is a haplotype pair for gir with i ≤ k. note that, the values of at wild card positions in r are still not yet determined. here we refer to hi,2r as a partially inferred haplotype segment on r. let intq={h r,h r,…,hk,2r}, where if two haplotype segments hi,2r and hi′,2r are identical, then we just keep one of them. note that if any haplotype segment hi,2r∈intq is undetermined at a position p, then all the haplotype segments in intq are undetermined at p. after partially determining hr and for every gir in dr, we use a heuristic method to infer the haplotype segments for gir∈nr. since we want to minimize the total number of resulting distinct haplotype segments on r, our strategy is to let the inferred haplotype segments for gir∈nr share as many haplotype segments as possible. this is actually clark’s inference rule. let q be a queue that contains a set of inferred haplotype segments on r. initially, q = intq. a partially inferred haplotype segment h in q can solvegjr if the following conditions hold: if h is at a position p then gjr is or at position p. if h is at a position p then gjr is or at position p. we can use h to solvegjr by constructing two haplotype segments as follows: using h to solve gjr: if h is at position p in r then we set hj,1r= at p and hj,2r at p is set according to rules c and c if h is at position p in r then we set hj,1r= at p and hj,2r at p is set according to rules c and c if h is undetermined at position p in r, and gjr is at p, then set hj,1r=h=hj,2r= at p. here we also have to determine the value of h at position p accordingly. after the undetermined value of h at p is determined, if h is obtained from a gir∈dr, then we also have to determine the values of hr and other hi,2r for gir∈dr at position p according to the haplotype value of h at p and rules c and c if h is undetermined at position p in r, and gjr is at p, then hj,1r and hj,2r remain undetermined at p. a genotype segment gir∈nr is solved if the pair of haplotype segments for gir are determined. in our algorithm, we use p to store the set of genotypes in nr that have not been solved. initially, p = nr. we then use the haplotype segments in q one by one and try to solve each of the genotypes in p. after trying to use a h ∈ q to solve all gjr’s in p, we delete h from q. two haplotype segments on r in q are compatible if there does not exist any position p in r such that one segment has value and the other segment has value at p. for two compatible haplotype segments h and h on r, we can merge them to form one haplotype segment h, where the value of h is determined as or if at least one of h and h is or and the value of h remains undetermined if both h1 and h are undetermined. again, for any position p, if the value of h1 or h is not identical to that of h, then the value of h1 or h is changed from undetermined to or thus, we have to update some of the previously inferred haplotype segments accordingly. when we use a h ∈ q to solve a gjr in p, we can obtain another new haplotype segment h′ on r. if h′ is compatible with a haplotype segment in q, we then merge them. note that, h′ might be compatible with more than one haplotype segment in q. in this case, we arbitrarily choose a compatible haplotype segment in q and merge the two haplotype segments. if h′ is not compatible with any haplotype segment in q, we add h′ into q. we give an example to illustrate the above process. example 1 dr = {g1 =  g2 =  g3 = 10221} and nr = {g4=10121}. after procedure p hr = 10∗11 and consequently intq = {h  = 10∗ h  = 11∗ h  = 10∗01}. after that, we can use h  = 10∗11 in intq to solve g4 = 10121 in nr. based on h  = 10∗ g4 = 10121 can be solved as h  = 10111 and h  =  moreover, since we want h and h to be identical, h is updated as h  =  correspondingly, we update hr = h  =  h  =  h  =  and h  =  after that, the set of distinct haplotype segments we have obtained so far is {hr =  h  = h  =  h  =  h  =  h  = 10101}. note that h and h are compatible, and the set of distinct haplotype segments is {hr =  h  = h  =  h  =  h  = h  = 10101}. after trying to use all the h’s in q to solve all gjr’s in p, q will become empty. when q is empty and p still contains at least two genotype segments, we consider all pairs of genotype segments gjr and gj′r in p and use the following method to infer the haplotype segments. inferring the haplotype segments from a pair gjr and gj′r in p: a position p in r is a conflicting position for gjr and gj′r if one of gjr and gj′r has the genotype value and the other has genotype value at p. the pair of gjr and gj′r can share a common haplotype segment on r if there is no conflicting position in r for gjr and gj′r. the shared haplotype segment can be computed as follows: if one of the genotype values at position p is then the haplotype value is at p; if one of the genotype values is at p, then the haplotype value is at p; if both genotype values are at p, then the haplotype value at p is undetermined. once the shared haplotype segment for gjr and gj′r are computed, we can determine the other haplotype segments for gjr and gj′r based on c and c after inferring the haplotype segments from a pair gjr and gj′r that can share a common haplotype segment, we delete gjr and gj′r from p, merge compatible inferred haplotype segments, and insert the newly obtained haplotype segments into q. once q is not empty, we can use haplotype segments in q to solve the genotype segments in p again. the process is repeated until p is empty. the detailed algorithm is given as algorithm algorithm 1 mutation region haplotype inference input: two sets of genotype segments dr={g1r,g2r,…,gkr} and nr={gk+1r,gk+2r,…,gnr} on r. output:true if there is a solution. type i false if two diseased haplotype segments are conflict at a position in r; type ii false otherwise. 1: compute the center haplotype hr as in procedure p 2: ifhr does not exist then 3: return type i false 4: else 5: set hi,1r=hr for gir∈dr; 6: compute hi,2r according to hi,1r and c and c for each gir∈dr; set q={h r,h r,…,hk,2r} and p = nr. 7: end if 8: whileq ≠ ∅ and p ≠ ∅do 9: delete a haplotype segment hi,2r from q; 10: ifhi,2r can solve gjrthen 11: use hi,2r to solve gjr. add the newly obtained haplotype segments into q and delete gjr from p. 12: end if 13: end while 14: if there are at least genotype segments in pthen 15: if there exists a pair of genotype segments gjr and gj′r that can share a haplotype segment on rthen 16: infer the haplotype segments of gjr and gj′r and insert them into q, goto line 17: end if 18: if any inferred haplotype segments for some gjr∈nr on r is identical to hrthen 19: return type ii false. 20: end if 21: fix hs,1r and hs,2r for each gsr∈p so that hs,1r≠hr and hs,2r≠hr. 22: if line line fails then 23: return type ii false 24: else 25: return true. 26: end if 27: end if the following is an example to illustrate the case when q becomes empty. example 2 dr = {g1 =  g2 = 11210} and nr = {g3 =  g4 = 11121}. after procedure p hr = 11110 and intq = { }. after trying all the h’s in intq to solve gi’s in nr, q becomes empty and p = nr  ={ }. in this case, we look at both and in p and infer a shared haplotype h  = h  = 11111 and the other two haplotype segments h  = 01111 and h  =  the algorithm for the whole chromosome for an input segment r on a chromosome, if algorithm returns true, then r is a valid segment. in order to get the mutation regions, we decompose the whole chromosome into a set of disjointed length l segments. for each segment, we run algorithm to test if the segment is valid. after finding all the valid segments, we repeatedly merge two valid segments into a long valid segment if the two segments are within 2l snps and algorithm returns type ii false on all the segments in the gap.after the above merging process, we obtain several long valid segments. for each such long valid segment [sb,se), we run algorithm on the three segments [sb− l,se + l), [sb− l,se + l) and [sb,se) and select the longest one ) which returns true. since we impose that algorithm returns type ii false for the segments in gaps in the merging process, we can always ensure that algorithm returns true for [sb,se). extending [b,e) to the left and right: after we obtain r = [b,e) as discussed above, we try to extend the segment [b,e) to the left and right. on the segment r=[b,e), we have inferred hi,1r and hi,2r for each gir∈dr∪nr. these hi,1r and hi,2r form a collection of disjointed sets h h …,hm, where each hk is a set of identical haplotypes in {hi,1r,hi,2r|i=,…,n} on r. we extend the segment [b,e) to the left and right by looking at each position p. we first try p = b − q for q = ,… and then p = e −  + q for q = ,…. for each hk, if there exist some hi,j in hk such that gi are and others are then every hi,j in hk should be. if there exists a hi,j in hk that has been set to the conflict value before, then we know that position p is a conflicting position and p should not be extended to be part of r and the extension process to the current direction should stop. otherwise, we set every hi,j in hk to be. if p is not a conflicting position, after setting hi,j in hk to be, we can determine the value of hi,j′ according to the value of hi,j at p, the value of gi at p and rules c1 and c again, we should test if such a value of hi,j′ is consistent with the value of hi,j′ determined before. if conflict exists, then p is a conflicting position and the extension process to the current direction should stop. let hi,j′∈hk′. if there is no conflict, we should also update the value of all h’s in hk′. this recursive process continues until no further change can be made. if all hi,j in hk are we set hi,j as undetermined. the extension process stops when we find conflicts in both directions. the extended region obtained from. after the extension process, our program reports all the mutation regions obtained in the algorithm. the complete algorithm to find the mutation regions on the whole chromosome is shown as algorithm algorithm 2 the algorithm for the whole chromosome input: two sets of genotype on the whole chromosome d = {g g …,gk} and n = {gk + gk + …,gn}. output: the detected mutation regions 1:decompose the whole chromosome into segments of length l =  2:for each segment do 3: use algorithm to test if the segment is valid. 4:end for 5:merge the valid segments to form longer segments. 6:for each segment [sb,se) obtained in line do 7: select the longest segment of [sb −  l,se + l), [sb −  l,se + l) and [sb,se) which will return true by calling algorithm denote it as [b,e). 8: extend [b,e) to get a candidate mutation region. 9:end for 10:output all the mutation regions obtained in the for loop of lines 6-  in this section, we first show some experiments on simulated data. we then give a real case study to show that our program can also handle real data. a discussion is given at the end of this section. experiments on simulated data in order to evaluate the performance of our method and the feasibility of the mathematical model proposed in this paper, we write a program in c++ to produce simulated data. the program takes a pedigree and the haplotype data for the whole chromosome of each founder in the pedigree as the input. it generates the haplotype data for the remaining individuals in the pedigree using the standard χ model for recombination with m equal to and according to the male/female averaged genetic map for chromosome downloaded from hapmap. also see. the haplotype data of a non-founder in the pedigree are generated to randomly inherit one strand of the four-strand chromatid bundle from each parent of the non-founder. a mutation point is selected uniformly at random from the snp sites of the chromosome. each diseased offspring is forced to inherit the strand with the mutation point and the normal offspring are forced to inherit the strand without the mutation point. in this way, we can guarantee that there is exactly one true mutation region. note that the true mutation region must be a shared mutation region. see implementation for the definition. moreover, since we know the haplotype data of all the individuals in the simulations, we can easily find the shared mutation regions. by definition, there may exist more than one shared mutation region. to generate the simulated data, we randomly chose some of the haplotype data for chromosome of unrelated japanese in tokyo and han chinese in beijing in the database of hapmap project as the haplotype data for each founder. there are about snp sites on this chromosome. recall that our program takes two sets of individuals d and n and their genotype data as input. after generating the haplotype data of each individual, we only use some of the individuals and their associated genotype data as the input of our program. to evaluate the performance of our method, we used different pedigrees to evaluate our algorithm. figures and are pedigrees of generations with to diseased individuals in the dashed rectangle at the bottom. those individuals in the dashed rectangle at the bottom of each pedigree are the input of our program. the correctly detected mutation regions are the intersection of the regions reported by the computer program and the true mutation region. here, precision is defined as the number of snps in the correctly detected mutation regions divided by the total number of snps in the regions output by the program. the value of recall is defined as the number of snps in the correctly detected mutation regions divided by the total number of snps in the true mutation region. so, if the value of recall is then all the snps in the true mutation region have been reported by the program. similarly, if precision is then all the reported snps are in the true mutation region. we performed experiments for each pedigree. since there are about snp sites on the chromosome, we set l =  for each region reported by our algorithm, we define a score as follows: let dh be the number of distinct haplotypes on this region and length =  the length of this region. then the score of this region is defined as score = ∗length, where n is the total number of input genotypes in d∪n. this score can balance the length of the mutation region and the number of distinct haplotype segments on the region. with the longer region and smaller dh, the score becomes higher. to illustrate the quality of our program, we report the results when our program reports the region with the highest score and the first three regions with the highest scores, respectively. in fact, our program does not need this score in the computation. the program simply reports all the mutation regions. see the genotype data error handling in discussion. the precision and recall on the experiments are shown in table only the genotype of the individuals in the latest generation of pedigree 1−4 are known in this experiment. several mutation regions may be detected by our algorithm. in table the results when our program reports the region with the highest score are shown in the columns under “one region”. the results when our program reports three regions with the highest scores are shown in the columns under “three regions”. the precision and recall are calculated based on the true mutation region, the reported region, and the intersection of the reported region and the true mutation region. the precision’ and recall’ are calculated by replacing the true mutation region with shared mutation regions. the column “time” indicates the average time of our program by running experiments on each pedigree. from table we can see that the values of recall are from % to % and the values of precision are from % to % in the four pedigrees. when the number of diseased individuals is increased, the values of both recall and precision are improved significantly. when there are or diseased individuals, the value of recall is more than 97%. that is, the program can report most of the snps in the true mutation region. in practice, one can often get the genotype data for the individuals of the latest two generations. thus, we study this case by looking at different input individuals based on the pedigrees in figures and now, we study different sets of input individuals in the latest two generations of pedigree these different sets of input individuals in the latest two generations in the pedigree are given in figure a square/circle with a slash indicates that such individual is not included as part of the input though the individual is used in generating the simulated data. for the rest of test, we performed experiments for each case and show the average values. table shows the results for the different sets of input individuals in figure the individuals in the latest and latest two generations are not distinguished in our algorithm. we just input the genotype for all the individuals without the slash. again, the setting is similar to that of table from table we can see that the values of recall for different inputs are close to 99% except for 2d-3fam- 2d-2fam- and 2d-2fam- where the input contains only or diseased individuals. comparing table and table we can see that more input individuals do help improve the values of precision and recall. we also performed similar experiments for pedigree 2-. the results are similar to that in table and are given in additional file we also tested the program using pedigrees containing and generations and diseased individuals, respectively, in the latest generation. the four pedigrees containing generations are shown in additional file 1: figure s figure s figure s and figure s in the additional file. the four pedigrees containing generations are shown in additional file 1: figure s figure s figure s and figure s in the additional file. again, the input individuals are the individuals in the dashed rectangle at the bottom of the pedigree. the experiment results for generations and generations are shown in table and table respectively. the settings of table and table are similar to that of table table and table show that the performance of our program for and generations is similar to that for generations. similar to the case of generations, for pedigrees with and generations, we also tested various cases when some individuals of the latest two generations are available as input individuals. the results for and generations are similar to that of generations. the detailed results are given in the additional file. a real case study to illustrate the usefulness of our program, we applied our method to a set of real data originally from the phase ii hapmap database and was studied in. in, the authors studied two ceu families ceph and ceph. they identified a segment on chromosome shared by the four individuals na na na and na for this set of data, there are totally snp sites on the chromosome after the unknown genotypes are eliminated from the database. there are totally snp sites between 107m and 110m on chromosome starting at the 122348-th snp site and ending at the 128866-th snp site. we applied our program with the six individuals in the two families ceph and ceph as input and set the four individuals na na na and na as diseased individuals. we set l =  our algorithm found several segments shared by the four diseased individuals. the lengths of all the reported segments are approximately snp sites except the longest ones. all these segments are shown in table table shows the starting and ending point of the segments, the number of distinct haplotypes on the segment, and the score for each segment. as there is only one conflicting position between segments and, we should consider such a conflicting position as a data error. therefore, segment should be the predicted mutation. this segment starts at the 124561-th snp site and ends at the 129451-th snp site. the details are shown in figure we can see that the shared segment found by plink in starts at the 122348-th snp site while the starting position of our reported segment is for the subsegment that we did not report, we found conflicting positions in this subsegment containing snp sites. however, on the segment containing snp sites reported by our program, there is only one conflicting position. this is strong evidence that the subsegment is not shared by all the four diseased individuals. we also looked at the segment with length of snp sites on the right of our reported segment, and found conflicting positions among the snp sites. we can see that the subsegments and on the left and right of our reported segment have approximately the same number of conflicting positions. discussion haplotype inference methods as discussed before, if the haplotype data for each input individual are known, the problem of finding the true mutation region is straightforward. currently, there are several population-based phasing methods that can give accurate haplotype segments. however, these methods can only phase a small number of snps effectively and take an extremely long time to infer the haplotype for the whole chromosome. lrp in can phase more than snps simultaneously within a reasonable time. however, it is still very slow for phasing snps of a whole chromosome. moreover, they cannot directly report the true mutation region for a set of input individuals. on the other hand, our program can complete the computation in less than seconds for about snps with about to individuals. related mutation region detection methods to our knowledge, all the existing software packages need a clearly given pedigree as part of the input. if the pedigree is not known, most of the software packages do not work. our algorithm deals with the case where the input individuals are closely related but the pedigree is not given. merlin is widely used for linkage analysis, where a pedigree is required as part of the input. it works well on snp data due to the use of sparse trees. however, it can only analyze pedigrees of moderate size. when the family size is big, a large memory space is needed and the computation cannot successfully be completed. as shown in, merlin cannot report the results for some pedigrees, e.g., p and p in, where there are less than input individuals. however, our program can deal with the cases where the number of input individuals is large. the first row in table of the additional file shows the results for input individuals. we can see that our method can give very high precision and recall in this case. therefore, our algorithm can handle some cases which cannot be handled by merlin. the rule-based algorithm in uses a set of heuristics for haplotype inference with a given pedigree. it can give very accurate results when the number of family members is large enough and for each nuclear family the genotype data for both parents are available. however, it does not work well when the genotype data of one of the parents are missing in the nuclear family. if the data for both parents are missing, it does not work. liden is an extended software package of the algorithm in. it focuses on handling the case where the genotype data for one of the parents in a nuclear family are missing for the entire chromosome. but it still does not work when the genotype data for both parents in a nuclear family are missing or the family pedigree is not given in the input. plink in and beagle in can identify the shared haplotype segment between two individuals based on population-based linkage analysis. but it cannot automatically identify the mutation region taking a set of individuals as the input, which is expected to give more precise prediction. the above real case study has illustrated this. genotype data error handling the real datasets often contain errors. handling the genotype data errors is an important issue in practice. for our program, we have a pre-process step to delete all the snps containing missing data. that is, if the genotype data for an input individual at an snp site are missing, we delete this snp site from the input. without this step, we cannot get reasonable results for the real case study. when the genotype data contain errors, it is hard to detect and correct them. the errors may affect our program’s results in two ways: an snp site in the true mutation region may become a conflicting position due to error; and the number of distinct haplotype segments to explain the genotype data is increased. when occurs, our score for the detected mutation regions becomes worse. when occurs once, our program reports two detected regions with the conflicting position in between. see and in table of the real case study. when this kind of error occurs many times, our program reports many regions separated by a few snps in the middle. when the user looks at the results of our program, it is possible to realize that the few snps between two closely located reported regions are due to errors. this is similar to other linkage analysis programs such as merlin, where each snp site has a score, and the user decides a region with high scores as the true mutation region. our program may work for the situation where the input individuals are from multiple families. our algorithm tries to find regions shared by all the diseased individuals. thus, as long as the diseased individuals from multiple families share the same haplotype segment on the true mutation region, our program should be able to find such region. even if the haplotype segments from different families on the true mutation region are slightly different, the program should be able to report several smaller regions with a few missing snps in the middle. again, it is possible for users to figure out the whole true mutation by ignoring the few missing snps in the middle.  we developed a software package for linkage analysis where the input individuals are closely related, but the pedigree is not known. we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed by dan gusfield. our heuristic algorithm simply uses clark’s inference rule to infer haplotype segments. experiments show that our program can give very high value of recall in various cases. this implies that the program can report more than 90% of the true mutation region. the value of precision varies from 29% to 90%. when the precision is 29%, the size of the reported regions is three times that of the true mutation region. this is still very useful for narrowing down the range of the disease gene. availability and requirements project name: mrd project homepage:http://www.cs.cityu.edu.hk/~wenjuacui/software/mutationregion/index.html. the source code is also available. operating system: platform independent programming language: java other requirements: java. or higher license: none any restrictions to use by non-academics: none competing interests the authors declare that they have no competing interests author’s contributions lw proposed the topic and ideas for algorithms, wc implemented the program, and both authors devised and developed the method and prepared the manuscript. both authors read and approved the final manuscript. supplementary material additional file 1 supplementary material: this file includes several figures and additional experimental results mentioned in the paper. it contains the different set of input individuals on pedigree 2- in the paper, the pedigrees containing and generations and, diseased individuals in the latest generation respectively. the tables show the results on the input of the above figures. click here for file acknowledgements this work is supported by a grant from the research grants council of the hong kong special administrative region, china. lusheng wang is the corresponding author.
the value of precision varies from 29% to 90%. this is still very useful for narrowing down the range of the disease gene location. to our knowledge, no algorithm can give good results for this case. we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed by dan gusfield. our heuristic algorithm uses clark’s inference rule to infer haplotype segments. linkage analysis is the first step in the search for a disease gene. in this paper, we study a very important case, where the sampled individuals are closely related, but the pedigree is not given. experiments show that our program performs well. linkage studies have facilitated the identification of several hundred human genes that can harbor mutations leading to a disease phenotype. the recall value is from 90% to 99% in various cases. to solve this problem, we first developed some heuristic algorithms for haplotype inference without any given pedigree. we ran our program both on the simulated data and a set of real data from the phase ii hapmap database. this implies that the program can report more than 90% of the true mutation regions. when the precision is 29%, the size of the reported regions is three times that of the true mutation region.
 linkage analysis is the first step in the search for a disease gene. linkage studies have facilitated the identification of several hundred human genes that can harbor mutations leading to a disease phenotype. in this paper, we study a very important case, where the sampled individuals are closely related, but the pedigree is not given. this situation happens very often when the individuals share a common ancestor or more generations ago. to our knowledge, no algorithm can give good results for this case.  to solve this problem, we first developed some heuristic algorithms for haplotype inference without any given pedigree. we propose a model using the parsimony principle that can be viewed as an extension of the model first proposed by dan gusfield. our heuristic algorithm uses clark’s inference rule to infer haplotype segments.  we ran our program both on the simulated data and a set of real data from the phase ii hapmap database. experiments show that our program performs well. the recall value is from 90% to 99% in various cases. this implies that the program can report more than 90% of the true mutation regions. the value of precision varies from 29% to 90%. when the precision is 29%, the size of the reported regions is three times that of the true mutation region. this is still very useful for narrowing down the range of the disease gene location. our program can complete the computation for all the tested cases, where there are about snps on a chromosome, within seconds. 