Identifying haplotype block structure by using ancestor-derived model and MDL principle
Hironori Fujisawa
1, Minoru Isomura
2, Shinto Eguchi
1, Masaru Ushijima
2, Satoshi Miyata
2, Yoshio Miki
2, and Masaaki Matsuura
21
Institute of Statistical Mathematics, Tokyo
2
Genome Center, Japanese Foundation for Cancer Research, Tokyo
This manuscript is a draft version.
ABSTRACT
A method for identifying a haplotype block structure is presented using an ancestor- derived model and MDL principle. The haplotype block structure is caused by the existence of ancestral haplotype and recombination event. Using this idea, the ancestor- derived model is constructed. The whole statistical model is proposed as a modifica- tion of the ancestor-derived model to furthermore treat some non-standard data. The haplotype block structure is identified using the optimal model selected by the MDL principle. The simulation study shows that the proposed method is powerful from the viewpoint of hotspot sensitivity and robust to mutation except near the edge of sequence. The proposed method was applied to two real data sets; the JFCR data and the 5q31 data of Daly et al. (2001). The analysis of the former data gave a clear explanation on the result of conventional association study of antitumor drug. The analysis of the latter data presented a similar haplotype block structure to in Daly et al. (2001).
Introduction
Variation in the human genome sequence plays a powerful role as a biological marker to
detect a disease-related gene. A typical example of the variation is the single nucleotide
polymorphism (SNP), which has been investigated in various regions of genome. Based
on SNP markers on chromosome 5q31, Daly et al. (2001) reported a block structure,
more precisely, that the SNP markers had a very strong correlation within each block
and presented a sudden decay of correlation between blocks. Such a structure is now
referred to as a haplotype block structure and the place of sudden decay is referred to
as a recombination hotspot. They furthermore pointed out that most haplotypes could
be interpreted as recombinants of only a few common haplotypes, which implies a low
haplotype diversity. Hereafter such common haplotypes are referred to as ancestral haplotypes. They also regarded the haplotype block structure as a new marker in place of SNP and then dramatically improved a traditional linkage analysis. Similar haplotype block structures have been discussed thereafter (Patil et al. 2001; Gabriel et al. 2002; Zhu et al. 2003; Kamatani et al. 2004).
An important problem is to make a rational and automatic method for identifying a haplotype block structure. The methods proposed in the past can be classified into four types. The first type aims at achieving a low haplotype diversity. Daly et al.
(2001) intended a locally minimal haplotype diversity based on haplotypic heterozy- gosity. Patil et al. (2001) focused on haplotype-tagging-SNPs. The calculation task was reduced through the dynamic programming (DP) algorithm suggested by Zhang et al. (2002). The second type aims at detecting a recombination hotspot. Gabriel et al. (2002) and Wang et al. (2002) used the coefficient of linkage equilibrium (LD) and the four-gamete test statistic as indexes of the degree of recombination, respectively.
The third type is an attempt to directly combine the above two ideas (Zhu et al. 2003;
Kamatani et al. 2004). The fourth type is based on a statistical model and minimum- description-length (MDL) principle (Rissanen 1978), which manages both haplotype diversity and recombination hotspot. The method proposed in this paper belongs to the fourth type.
Koivisto et al. (2003) employed a simple statistical model, which is determined only by the haplotype block structure. The MDL principle selected the optimal model, which identifies the haplotype block structure. Greenspan and Geiger (2003) used a hidden markov model with the ancestral haplotype and mutation as well as the haplotype block structure. Their model is more faithful to a genetical property, but too complicated to be optimized. Anderson and Novembre (2003) discussed a markov model with the haplotype block structure, but without the ancestral haplotype and mutation. They incorporated a special probabilistic structure related to recombination event into the transition of the markov model, which was one of the keys of their paper.
They put computability ahead of faithful modeling, as is well-known that in this case the calculation task can be reduced through the DP algorithm. They showed that their method worked remarkably well in comparison with past methods in simulation study and real data analysis. Their method is referred to as the AN method in this paper.
This paper constructs a statistical model as follows. The most important reason
of generating the haplotype block structure is the existence of ancestral haplotype
and the recombination event. We first incorporate a probabilistic structure of recom-
bination event with ancestral haplotype into the central part of the model, which is
referred to as an ancestor-derived model, but the resulting model can treat only recom-
binants of ancestral haplotypes. A rest of the problem is the modeling of the remaining
non-recombinant haplotypes. They are not always interpreted by a simple genetical
property (e.g. mutation) and their frequency is small in general, so that we adopt a
simple full frequency model with no structure to avoid extra complexity of the model.
The objective statistical model is proposed as a mixture of ancestor-derived model and simple full frequency model. The parameter of the mixture model is estimated by the maximum likelihood principle through the EM algorithm (McLachlan and Krishnan 1997). The DP algorithm can be applied to reduce the calculation task.
The statistical model and parameter estimation have been discussed already. The rest part of the method for identifying the haplotype block structure is the selec- tion of the optimal model, so that the optimal model determines the haplotype block structure. This paper adopts the MDL principle as the model selection criterion. As described later in detail, the MDL principle prefers a larger maximum log-likelihood and a smaller number of underlying haplotypes, more abstractly, a better probabilistic structure fitting to the data and a lower haplotype diversity. This is much suitable for the purpose of identifying the haplotype block structure. The proposed method was implemented in the program ADBlock, which is available at the web site of ADBlock.
The ancestor-derived model is the same as a certain restricted markov model and can interpret recombinant haplotypes by a smaller number of parameters. This fact implies that the proposed method tends to be more powerful from the viewpoint of hotspot sensitivity than the AN method. This was verified in simulation study. The robustness to mutation was also investigated.
The proposed method was applied to two real data sets. One is the JFCR data, which was sampled in the Japanese Foundation for Cancer Research (JFCR), based on informed consent. The proposed method and the AN method were applied to the data and then the haplotype block structure was identified. The proposed method gave a clearer haplotype block structure than the AN method in relation to a prior association study to an adverse effect of antitumor drug. Since the number of individuals was only 75, the higher power of the proposed method would be useful. Another is the 5q31 data analyzed by Daly et al. (2001). The proposed method was applied to the data and the resulting haplotype block structure was similar to the conventional structure.
Methods
Ancestor-derived model
Some notations are prepared at first. Let the numbers of SNP markers and observed haplotypes be denoted by L and N . The i-th observed haplotype can be expressed as h
i= (h
i1, . . . , h
iL) for i = 1, . . . , N , where h
il∈ { 0, 1 } because the SNP is biallelic.
The suffix i is sometimes omitted for simplicity. The block structure B = { B
1, . . . , B
K} can be expressed as the partition of SNP numbers, where K is the number of blocks, B = { 1, . . . , L } = B
1∪ · · · ∪ B
K, B
k= { l
k−1+ 1, · · · , l
k} is the set of adjacent SNP numbers with l
0= 0 and l
K= L. The haplotype corresponding to the above partition is denoted by h = ³ h
(1), . . . , h
(K)´ , where h
(k)= ³ h
lk−1+1, · · · , h
lk´ is the partial haplotype
on block B
k. Let G = { g
1, . . . , g
A} be the set of ancestral haplotypes and q
athe
frequent probability of g
a, where P
Aa=1q
a= 1. Assume that the observed haplotype h is a recombinant haplotype, more precisely, h = ³ g
a(1)1, . . . , g
(K)aK´ for some (a
1, . . . , a
K).
Let R = (R
1, . . . , R
K−1) be the latent indicator whether the recombination event happens or not, where R
k= 1 if the recombination event happens between two blocks B
kand B
k+1and R
k= 0 otherwise. Let the recombination rate be denoted by λ
k= Pr(R
k= 1). Assume that R
k’s are independent variables. The probability of R is expressed as
Pr(R = r) =
K
Y
−1k=1
λ
rkk(1 − λ
k)
1−rk.
Suppose that R is given in the following. Another partition of the SNP numbers is denoted by B = B
[1]∪ · · · ∪ B
[K∗], where K
∗= P
Kk=1−1R
k+ 1, B
[k∗]is the union of some adjacent B
(k)’s, the recombination events happen only between B
[k∗]’s and not within each B
[k∗]. A relation between two block partterns is displayed in Figure 1. The haplotype corresponding to the above partition is denoted by h = ³ h
[1], . . . , h
[K∗]´ . The partial haplotypes h
[k∗]’s derive from g
a’s because the haplotype is a recombinant of ancestral haplotypes. This correspondence can be expressed as the indicator variable C = (C
[1], . . . , C
[K∗]), where C
[k∗]= a means that h
[k∗]derives from g
a. Note that the variable C depends on the structure of R. The conditional probability of C given R is defined as
Pr(C = c | R = r) =
K∗
Y
k∗=1
Y
Aa=1
q
aI(C[k∗]=a),
where I( A ) is one if A is true and zero otherwise. Consequently, the probability that the latent variable is observed is given by Pr(R = r, C = c) = Pr(R = r) Pr(C = c | R = r), which is called a complete model (McLachlan and Peel 2000).
Recall the original problem. The event where the haplotype h is observed is ex- pressed by a latent variable as F
h= n (R, C) | h = ³ g
C[1][1], . . . , g
C[K[K∗]∗]´o , which implies the objective frequent probability
Pr(H = h) = X
(r,c)∈Fh
Pr(R = r, C = c),
which is referred to as the ancestor-derived model. This model is determined by the haplotype block structure, B , and the set of ancestral haplotypes, G . The parameter of the model consists of the recombination rate, λ = (λ
1, . . . , λ
K−1), and the frequency of ancestral haplotype, q = (q
1, . . . , q
A).
To concretely understand the ancestor-derived model, illustrate the simple case
where the number of SNP markers is five, L = 5, the block structure is given by
B
1= { 1, 2, 3 } and B
2= { 4, 5 } , the number of ancestral haplotypes is two, A = 2,
given by g
1= (0, 0, 0, 0, 0) and g
2= (1, 1, 1, 1, 1). Suppose that the observed haplotype is h = (0, 0, 0, 1, 1). The recombination event certainly happens between two blocks, R = 1. The observed haplotype is a result of combination of two partial descendents, more precisely, h = (h
[1], h
[2]) = (g
1[1], g
2[2]) and C = (1, 2). It therefore follows that F
h= { (1, (1, 2)) } and
Pr(H = h) = Pr(R = 1, C = (1, 2)) = Pr(R = 1) Pr(C = (1, 2) | R = 1) = λq
1q
2. Suppose that the observed haplotype is h = (0, 0, 0, 0, 0). If no recombination event happens, R = 0, then only one block is present, K
∗= 1, and the observed haplotype is a direct descendent of the ancestral haplotype g
1, C = 1. Consider the case where the recombination event happens, R = 1. The number of blocks is two, K
∗= 2.
The observed haplotype is a result of combination of two partial descendents, more precisely, h = (h
[1], h
[2]) = (g
1[1], g
1[2]) and C = (1, 1). It therefore follows that F
h= { (0, 1), (1, (1, 1)) } and
Pr(H = h) = Pr(R = 0, C = 1) + Pr(R = 1, C = (1, 1))
= (1 − λ)q
1+ λq
12. A general case can be extended by a similar way.
Mixture model and parameter estimation
The ancestor-derived model has been constructed to treat recombinant haplotypes.
A rest of the problem is the modeling of the remaining non-recombinant haplotypes.
As described already in Introduction, a simple full frequency model is applied to their haplotypes. Let U = { u
1, . . . , u
D} be the set of distinct non-recombinant haplotypes and p
dthe frequent probability of u
d, where P
Dd=1p
d= 1. The full frequency model is given by Pr(H = h) = Q
Dd=1p
I(h=ud d). The objective model is proposed as a mixture of ancestor-derived model and simple full frequency model with a mixing proportion ω.
The parameter of the mixture model, θ = (λ, q, p, ω), can be estimated by the maximum likelihood principle. Let the maximum likelihood estimate of θ be denoted by θ. Note that the observed haplotype belongs to either of two underlying models b and never to both simultaneously. Let the set of recombinant haplotypes be denoted by { h
1, . . . , h
n} . It is clear that ˆ ω = n/N , ˆ p is simply given by the observed frequency, and ˆ ξ = (ˆ λ, q) is the maximizer of ˆ P
ni=1log Pr(H = h
i; ξ) with respect to ξ = (λ, q), which can be obtained through the EM algorithm (McLachlan and Krishnan 1997).
The exact expression is derived in Appendix A.
When the frequent probability, Pr(H = h
i), is evaluated, much calculation task
may be necessary, because the size of F
his an exponential order of the number of
blocks, more precisely, (A + 1)
K−1. This task can be reduced up to the polynomial
order of low degree by virtue of the DP algorithm. A detailed definition of the algo-
rithm is given in Appendix A. In particular, its effectiveness appears in the parameter
estimation because the frequent probability is calculated for each iteration step of the EM algorithm.
Code length of the model
Remember that the mixture model is determined by the block structure, B , the set of ancestral haplotypes, G , the set of the distinct non-recombinant haplotypes, U , and the resulting probability structure. By summing up their code lengths, we can obtain the whole code length of the model. In the following, the code length is concretely calculated. The optimal model is selected as the minimizer of the code length by the MDL principle. For detailed discussion, see e.g. Hansen and Yu (2001).
The block structure can be determined by the number of blocks, K, and the cor- responding block patterns. The former ranges from 1 to L. Let ν
L,K= (L − 1)!/(K − 1)!(L − K)!, which is the number of possible block patterns. Hence, the necessary length of coding the haplotype block structure is given by log L + log ν
L,K. Next consider two sets G and U , which are determined by their sizes and the sequences of haplotypes.
The size of the combined set ranges from 1 to 2
L, whose necessary code length is log 2
L= L. From this code length, the size of the combined set is known as A + D. We must furthermore know the size of the first set, A, to determine each size of two sets.
Because the maximum number of this size is already known as A + D, the necessary length of coding this size is given by log(A + D). The code length L is necessary to express a sequence of haplotype because each component is biallelic and the number of components is L. The necessary length of coding all the haplotypes of G and U is given by (A + D)L, because the number of haplotypes is A + D. Consequently, the necessary length of coding two sets is given by (A + D + 1)L + log(A + D). The code length of the probability structure is approximated by minus maximized log-likelihood plus log N/2 times the number of parameters. The total number of parameters is (A + D − 2) + (K − 1) + 1 = A + D + K − 2. It therefore follows that the whole code length of the model is given by
ψ = log L + log ν
L,K+ (A + D + 1)L + log(A + D)
−
X
Ni=1
log Pr(H = h
i; θ) + (A b + D + K − 2) log N/2.
There exist two problems when the MDL principle is directly applied to the real
data. The existence of the very rare haplotype may lead to an erroneous conclusion,
which is well-known as over-sensitivity to outlier in a statistical field. The number of
distinct observed haplotypes is many in general and hence the number of candidates
of G is extraordinarily large. The former problem can be avoided by not using the
very rare haplotype. The latter problem can be overcomed using a prior information
that the ancestral haplotype has a very larger frequency. For a detailed procedure, see
Appendix B.
Comparison of the ancestor-derived model with the markov model
The ancestor-derived model can interpret aecombinant haplotypes by a smaller number of parameters than the markov model constructed by Anderson and Nobe- mbre (2003). The markov model has a flexible transition structure with A strength parameters to connect A ancestral haplotypes at each recombination hotspot. The ancestor-derived model has only one parameter as a recombination rate at each recom- bination hotspot. It follows that the extra number is (A − 1)(K − 1), which implies that the proposed method has a higher possibility of identifying a correct haplotype block structure. This was verified in simulation study.
More than one generation change
Let Pr
[m](H = h) be the probability that the haplotype h is observed at m-th generation as a recombinant of the ancestor haplotypes of 0-th generation. The case m = 1 corresponds to the ancestor-derived model. For simplicity, consider the case K = 2. Let the recombination rate at each generation change be denoted by ζ. It is shown that
Pr
[m](H = h) = (1 − ζ)
mPr
[0](H = h)
+ { 1 − (1 − ζ)
m} Pr
[0]³ H
(1)= h
(1)´ Pr
[0]³ H
(2)= h
(2)´ . The derivation is given in Appendix C. It is easily verified that this is completely the same as the ancestor-derived model where the recombination rate λ is replaced by { 1 − (1 − ζ)
m} . Thus, the ancestor-derived model has a reasonable understanding.
Furthermore, since ζ is sufficiently small in general, it holds that λ ≈ mζ, so that the recombination rate λ of the ancestor-derived model approximately corresponds to the true recombination rate times the number of generation changes.
Results
Simulation Study
An artificial sample was generated as follows. Two underlying sets of ancestral
haplotypes and their frequencies are given in Table 1, where the case A = 4 is the
same as the blocks 4 and 5 of Daly et al. (2001). The recombination hotspot was set
between SNP markers 5 and 6 with the recombination rate λ = 0.1, 0.3. Each setting
uniquely determines the ancestor-derived model. We randomly sampled n haplotypes
from the model for n = 200, 500, 1000 and then exposed each component of haplotype to
mutation with rate µ = 0.01, 0.005, 0.001. Note that their mutation rates correspond to
the cases where the rate of non-exposed haplotype is (1 − µ)
10≈ 1 − 10µ = 0.9, 0.95, 0.99,
respectively. We drew upon Daly et al. (2001) to set various parameter values.
Based on 100 samples, the proposed method and the AN method were compared from the viewpoints of the hotspot and non-hotspot sensitivities, which are the rates of correct judgment of hotspot and non-hotspot. The case where the estimated recom- bination rate was less than 0.03 was neglected because such a case does not present a sufficient evidence of recombination hotspot. The results of the case µ = 0.01 are displayed in Tables 2 and 3, where the spot k means the place between SNP markers k and k + 1. The hotspot and non-hotspot sensitivities corresponds to the spot 5 and the other spots, respectively. The case µ = 0.001 presented an almost full efficiency from the viewpoint of non-hotspot sensitivity and the case µ = 0.005 was intermediate (not shown).
The proposed method was superior to the AN method from the viewpoint of hotspot sensitivity, especially in the case where n = 200. The reason has been described already.
Focus on the case where A = 2 and n = 1000. The proposed method was more robust to mutation than the AN method. The reason will be that the markov model was too flexible. On the other hand, the proposed method was less stable than the AN method from the viewpoint of non-hotspot sensitivity near the edge of sequence.
A typical example of incorrectly regarding non-hotspot as hotspot is the following.
Suppose that the haplotype data include two ancestral haplotypes 00 · · · 0 and 11 · · · 1 and mutant haplotypes (not recombinant) 10 · · · 0 and/or 01 · · · 1. Forcedly assume that the first spot is a hotspot. The maximum log-likelihood will decrease because of erroneous modeling, but the degree of decrease may not be so large because the SNP marker is at most biallelic. On the other hand, the necessary length of coding U decreases because the mutant haplotypes can be expressed as recombinants of two ancestral haplotypes and then the mutant haplotypes are taken away from U . Such a trade-off of code length could arises a decrease of the whole code length, which led to an incorrect haplotype block structure. A suspicious hotspot near the true hotspot was also explained by the same reason.
As a conclusion, the proposed method is powerful from the viewpoint of hotspot sensitivity and stable near the center of sequence, but not always powerful near the edge of sequence.
JFCR data
The JFCR investigates associations between SNP markers and adverse effects of
antitumor drugs. Recently, an association between some SNP markers and a certain
adverse effect has been detected using the standard association study. The layout of
SNP markers is displayed in Figure 2 and the resulting p-values are given in Table 4. (A
detailed result will be published elsewhere.) It was considered that the SNP markers
7-9 were objective markers because the p-values were sufficiently small and the SNP
markers 3-6 were relevant markers because the p-values were relatively small. The
sample size was 75. The proposed method was applied to the JFCR data and then the
structure identified was much useful to interpret the result of association study.
The haplotype data were restored from the genotype data by the Haplotyper (Niu et al. 2002), which is based on the model-based approach (Excofficer and Slatkin 1995) and the Bayesian inference. Using the restored haplotype data, the proposed method gave the haplotype block structure B = {{ 1, 2 } , { 3 } , { 4, 5, 6 } , { 7, 8, 9 } , { 10 }} with re- combination rates 0.145, 0.474, 0.047, and 0.160. The number of ancestral haplotypes was four, but on the fourth block the number of common haplotypes was only two. It was expected that the fourth block is a significant region associated with the adverse ef- fect and that the third and fourth blocks are strongly related because the recombination rate between their two blocks was small. A further research has been done. The AN method identified the haplotype block structure B = {{ 1, 2, 3 } , { 4, 5, 6, 7, 8, 9 } , { 10 }} , which missed the spots 2 and 6. The missing of spot 6 may be due to a lack of power for a small sample size.
5q31 data
Daly et al. (2001) reported the haplotype block structure on chromosome 5q31 with 103 SNP markers. Using 258 genotypes, Anderson and Novembre (2003) identified the haplotype block structure similar to in Daly et al. (2001). Using the same data, the proposed method also identified the similar structure.
The genotype data include many missing responses. It may be possible to fill the missing responses by some alternatives, but uncertainty of the filled responses leads to an unreliable conclusion. For this reason, the genotype whose missing rate is more than 10% was not used. The haplotype data were restored from the genotype data by the software Haplotyper. The proposed method is powerful near the center, as seen in the simulation study, and needs much calculation task when the number of SNP markers is large. Therefore, the following procedure was adopted to identify the haplotype block structure. First the proposed method was applied to the first 16 SNP markers 1-16, including 15 spots, and we determined whether the first 10 spots 1-10 are hotspots or not. Next the proposed method was applied to the next 16 SNP markers 6-21 and we determined whether the next 5 spots 11-15 are hotspots or not. This step was continued and then the last 7 spots were determined using the remaining SNP markers 91-103.
Focus on the first 16 SNP markers to illustrate how the proposed method worked.
The haplotype block structures with smaller code lengths are given in Table 5. The smallest code length shows that the only one recombination hotspot exists at the spot 8. The second case implies a possibility that the spot 9 is a hotspot in place of the spot 8, because the second code length is not so far from the smallest code length.
The third to sixth cases imply two possibilities. One is the existence of hotspot at the
spot 14 or 15 and the other is that the proposed method is not stable near the edge of
sequence. The seventh case presents suspicious adjacent hotspots, as described in the
simulation study. The last case corresponds to the one where no hotspot is present, although the corresponding code length is much far from the smallest code length. The identified ancestral haplotypes and their estimated frequencies are displayed in Table 6.
On the left side of the hotspot, only two haplotypes are observed. On the right side, the haplotypes are similar except for the third. The fifth haplotype may appear as a mutant of the second haplotype at a very old era and thereafter establish itself. The spot 15 was identified as a hotspot with the large recombination rate λ = 0.3 by the next analysis and hence the fourth haplotype may be regarded as a recombinant.
Two sets of hotspots identified by the proposed method and Daly et al. (2001) were given by
{ 8,15,28,37,39,44,86,90,91,92,98 } ,
{ 8/9,14/15,24,35,40,45,76/77,84/85,91,98 } ,
where e.g. 8/9 may mean that either/both is/are the hotspot (because they did not clearly identify the hotspots in their paper). The two sets are similar except for two different points. One is the missing of the hotspot 76/77 and the other is the existence of adjacent hotspots from the spot 90 to 92. The missing of the hotspot 76/77 may be caused by the difference of the methods for restoring genotype data to haplotype data. In our restored haplotype data, the number of recombinant events illustrated in Daly et al. (2001) was only one. The reason of generating adjacent hotspots was the same as that described in the simulation study. The spot 91 will be the true hotspot and the spots 90 and 92 will be suspicious hotspots. The reason is as follows: The recombination rates of the spots 90, 91, and 92 were estimated as 0.041, 0.271, and 0.033, respectively. The first and third recombination rates were not so large. The spot 91 was always identified to be a hotspot among the haplotype block structures with the top 20 code lengths, but the spots 90 and 92 were not.
Discussion
This paper focuses on the existence of ancestral haplotype and the probabilistic structure of recombination event to construct the statistical model. The haplotype block structure is identified using the optimal model selected by the MDL principle.
The simulation study shows that the proposed method is powerful from the viewpoint of hotspot sensitivity and robust to mutation except near the edge of sequence. Two real data analyses implies that the proposed method works well. In the following, some issues for the future are discussed.
Some tuning parameters are prepared when the proposed method is applied, as
described already (see Appendix B). To avoid over-sensitivity of the model, the hap-
lotype whose observed number was one or two was not used in data analysis. When
such a device was not done, we often identified suspicious hotspots. To reduce the
caluculation task, the set of ancestral haplotypes was constructed using two threshold values; 0.1 and 0.03. When the threshold values were slightly changed to 0.2 and 0.01, a major change of the result was not observed. In the simulation study, the small recombination rate was neglected. In a real data analysis, it is recommended that the small recombination rate should be reviewed comparing with the original data set. The user of the software ADBlock can change tuning parameters, depending on data.
The mutation event is an important factor to understand genome sequences but is not incorporated into the mixture model to avoid extra complexicity of the model. In fact, incorporating the mutant event into the model is worth being studied, but it would be much difficult to optimize this model because this model causes heavy calculation task and over-sensitivity to small frequencies relevant to the mutation event.
Remember that the haplotype and the genotype has a latent relation. We can prepare a new latent variable to introduce its latent relation in addition to the latent variable (R, C ). This new latent variable will enable us to construct an extended model to directly treat the genotype data. To make this idea feasible, it will be necessary to overcome another heavy task problem.
The proposed method identifies the optimal haplotype block structure but presents no statistical significance of the optimal structure against other candidates. A simple and feasible method may be to calculate the confidential probability of the optimal model, based on the bootstrap sampling (Hall 1992). This is an important issue and a further detailed research will be desirable.
Appendix A
Let us give some notations before derivations. Let h
(k1;k2)and h
[k∗1;k∗2]be the haplo- types on B
(k1)∪ · · · ∪ B
(k2)and on B
[k1∗]∪ · · · ∪ B
[k2∗], respectively. Let R
(k)and R
[k∗]be the vector of the recombination event on blocks B
(k)and B
[k∗]and their boundaries, respectively, and let R
(k1;k2)and R
[k1∗;k2∗]be the vector whose components are the recom- bination variables within B
(k1)∪ · · · ∪ B
(k2)and B
[k∗1]∪ · · · ∪ B
[k2∗]except for boundaries, respectively.
EM algorithm
The EM algorithm is expressed as the following iteration step from ˆ ξ
(t)to ˆ ξ
(t+1): ξ ˆ
(t+1)= arg max
ξ
E
(t)"
nX
i=1
log Pr(C = C
i, R = R
i; ξ) | H = h
i#
,
where the expectation E
(t)[ · ] means E h · ; ˆ ξ
(t)i . The convergence value of ˆ ξ
(t)is the maximum likelihood estimate. The logarithm of the complete model is given by
log Pr(C = c, R = r) = log Pr(R = r) + log Pr(C = c | R = r)
=
K
X
−1k=1
{ r
klog λ
k+ (1 − r
k) log(1 − λ
k) } +
X
Aa=1 K∗
X
k∗=1
I(C
[k∗]= a) log q
a. Let the conditional expectations be denoted by
r
ki(t)= E
(t)[r
k| H = h
i] , s
(t)ai= E
(t)"
K∗X
k∗=1
I(C
[k∗]= a) | H = h
i#
,
and let the sample means be denoted by ¯ r
k(t)= P
ni=1r
ki(t)/n and ¯ s
(t)a= P
ni=1s
(t)ai/n. The iteration algorithm is re-expressed as
ξ ˆ
(t+1)= arg max
ξ
X
Kk=1
n r ¯
k(t)log λ
k+ (1 − r ¯
(t)k) log(1 − λ
k) o +
X
Aa=1
¯
s
(t)alog q
a, which implies that ˆ λ
(t+1)k= ¯ r
k(t)and ˆ q
a= ¯ s
(t)a.
Let Pr
(t)( · ) = Pr ³ · ; ˆ ξ
(t)´ . It holds that
Pr(H = h | R
k= 1) = Pr(H
(1;k)= h
(1;k)| R
k= 1) Pr(H
(k+1;K)= h
(k+1;K)| R
k= 1)
= Pr(H
(1;k)= h
(1;k)) Pr(H
(k+1;K)= h
(k+1;K)),
because under R
k= 1 the left-hand and right-hand sides are independent by virtue of the latent variable model of C given R, and the probability of the observed haplotype is independent of the recombination event outside of the observed part. Similar identities are often used in the following without any special description. It follows that
r
ki= E
(t)[r
k| H = h
i]
= Pr
(t)(r
k= 1 | H = h
i)
= Pr
(t)(r
k= 1) Pr
(t)(H = h
i| r
k= 1) . Pr
(t)(H = h
i)
= ˆ λ
(t)kPr
(t)³ H
(1;k)= h
(1;k)i´ Pr
(t)³ H
(k+1;K)= h
(k+1;K)i´ . Pr
(t)(H = h
i) and that
s
(t)ai= E
(t)"
K∗X
k∗=1
I(C
[k∗]= a) | H = h
i#
= X
r K∗
X
k∗=1
Pr
(t)³ C
[k∗]= a, R = r | H = h
i´
= X
r K∗
X
k∗=1
Pr
(t)³ H = h
i| C
[k∗]= a, R = r ´ Pr
(t)³ C
[k∗]= a | R = r ´
Pr
(t)(R = r) . Pr
(t)(H = h
i)
= X
r K∗
X
k∗=1
Pr
(t)³ H
[1;k∗−1]= h
[1;ki ∗−1]| R
[1;k∗−1]= r
[1;k∗−1]´ Pr
(t)³ H
[k∗]= h
[ki ∗]| C
[k∗]= a, R
[k∗]= (1, 0, . . . , 0, 1) ´ Pr
(t)³ H
[k∗+1;K∗]= h
[ki ∗+1;K∗]| R
[k∗+1;K∗]= r
[k∗+1;K∗]´ ˆ
q
a(t)Pr
(t)³ R
[1;k∗−1]= r
[1;k∗−1]´ Pr
(t)³ R
[k∗]= r
[k∗]´ Pr
(t)³ R
[k∗+1;K∗]= r
[k∗+1;K∗]´ . Pr
(t)(H = h
i)
= X
r K∗
X
k∗=1
Pr
(t)³ H
[1;k∗−1]= h
[1;ki ∗−1]R
[1;k∗−1]= r
[1;k∗−1]´ Pr
(t)³ H
[k∗+1;K∗]= h
[ki ∗+1;K∗]R
[k∗+1;K∗]= r
[k∗+1;K∗]´ I ³ h
[ki ∗]= g
a[k∗]´ q ˆ
a(t)Pr
(t)³ R
[k∗]= r
[k∗]´ . Pr
(t)(H = h
i)
= ˆ q
a(t)X
k1<k2
Pr
(t)³ H
(1;k1)= h
(1;ki 1)´ Pr
(t)³ H
(k2+1;K)= h
(ki 2+1;K)´ I ³ h
(ki 1+1;k2)= g
a(k1+1;k2)´
Pr
(t)³ R
k1= 1, R
(k1+1;k2)= (0, . . . , 0), R
k2= 1 ´ . Pr
(t)(H = h
i)
= ˆ q
a(t)X
k1<k2
Pr
(t)³ H
(1;k1)= h
(1;ki 1)´ Pr
(t)³ H
(k2+1;K)= h
(ki 2+1;K)´ I ³ h
(ki 1+1;k2)= g
a(k1+1;k2)´ λ ˆ
(t)k1k
Y
2−1k0=k1+1
³
1 − λ ˆ
(t)k0´ λ ˆ
(t)k2. Pr
(t)(H = h
i) .
DP algorithm
Let v
k= ( · , . . . , · , 1, 0, . . . , 0) be the (K − 1)-dimensional vector, where the k-th component is one, the components of the right-hand side are zero, and those of the left-hand side are zero or one. Note that any variable of the recombination event can be expressed as a certain v
kand that v
k(1;k)is completely free and v
k(k;K)= (1, 0, . . . , 0).
The frequent probability of the observed haplotype is expressed as Pr (H = h) =
K
X
−1k=0
Pr (H = h, R = v
k)
=
K
X
−1k=0
Pr (H = h | R = v
k) Pr (R = v
k)
=
K
X
−1k=0
Pr ³ H
(1;k)= h
(1;k)| R
(1;k)= v
k(1;k)´
Pr ³ H
(k+1;K)= h
(k+1;K)| R
(k;K)= v
k(k;K)´ Pr ³ R
(1;k)= v
k(1;k)´ Pr ³ R
(k;K)= v
k(k;K)´
=
K
X
−1k=0
Pr ³ H
(1;k)= h
(1;k)´
X
Aa=1
q
aI ³ g
(k+1;K)a= h
(k+1;K)´ λ
kK
Y
−1k0=k+1
