Statistical models for studying DNA sequence evolution

PHYSICA/ ELSEVIER

Physica A 221 (1995) 159-167

Statistical models for studying DNA sequence evolution Wen-Hsiung

Li, Xun Gu

Human Genetics Center, SPH, University of Texas, PO. Box 20334, Houston, TX 77225, USA

Abstract Statistical models for DNA sequence evolution and several methods for estimating the number of substitutions per site between two sequences are discussed. These methods are modified to include the effect of rate variation among sites.

1. Introduction The standard distance measure between two DNA sequences has always been the number of substitutions per site ( K ) , for several reasons. First, K is often used to reconstruct phylogenetic trees (i.e. the evolutionary relationships among extant species) [ 1 ]. Second, K is very useful for studying the pattern and mechanism of DNA evolution, for example, for testing the molecular clock hypothesis [2] and for detecting selection in sequence evolution [3]. Third, under some assumptions, K can be used to estimate the divergent time (T) between species [4]. However, K is generally not equal to the observed number of differences per site between two DNA sequences, because multiple substitutions at a given site may occur, especially when the sequence divergence is large. Therefore, to estimate K, a stochastic (Markov) model for DNA evolution is required. In this paper, we discuss such stochastic models and methods used to estimate K.

2. Models of nucleotide substitution Since at each site in a D N A sequence, there are four possible nucleotides (A, T, C and G), a model of nucleotide substitution is characterized by a 4 × 4 rate matrix R, which is also called the pattern of nucleotide substitution. This matrix can be represented as follows: 0378-4371/95/$09.50 (~) 1995 Elsevier Science B.V. All rights reserved SSD! 0378-437 1 (95) 00235-9

160

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167 A A T C G

T

-(a+b+c) d g j

C

a -(d+e+f) h k

G

b e -(g+h+i) l

c f i -(j+k+l)

The ijth element of R, rij, is the substitution rate from nucleotide i to nucleotide j, i 4: j; the diagonal elements are given by rii = - ~ j , i rij, so that the sum of the elements in each row is zero. Thus, the most general model has 12 independent parameters but this model is too complex to apply. Therefore, it is necessary to make some assumptions on R to develop a useful method for estimating K. As examples, the following are two simple but widely-used models. The first one is Jukes and Cantor's one-parameter model [5], in which the substitution rate is assumed to be the same for all nucleotides, i.e. r~i = u, for all i 4: j. For this model, the rate matrix R is

A T C G

A

T

C

G

-3u u u u

u -3u u u

u u -3u u

u u u -3u

The second is Kimura's two-parameter model [6], in which the rate of transitional substitution (i.e. changes between A and G or between C and T) may not be equal to that of transversional substitution (i.e. all the other types of changes) ; the two rates are denoted by s and v, respectively. This model is more realistic than the one-parameter model because transitions are generally more frequent than transversions in DNA evolution. The rate matrix R under this model is

A T C G

A

T

C

G

- ( 2 v + s) v v s

v - ( 2 v + s) s v

v s

s v v

- ( 2 v + s) v

- ( 2 v + s)

In general, let Pij(t) be the transition probability from nucleotide i to j after t time units, e.g., for nucleotide A at time 0, PnG(t) is the probability of having G after time t. By Markovian theory, Pij(t) satisfies the linear differential equation

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167

×

161

Y

Fig. 1. Two sequences diverged t time units ago.

dPij (t)

-~-~rikPkj(t),

-fit

i,j = A,G,T,C

(1)

k

with the initial condition P/i(0) = 1 and P/j = 0 (i 4: j ) . Let matrix P ( t ) consist of Pij(t). In matrix form, (1) can be expressed as dP(t) dt

RP(t)

-

(2)

with the initial condition P ( 0 ) = I, where I is the identity matrix. The solution of (2) is P ( t ) = e Rt.

(3)

Because we generally do not know the initial sequences, to study DNA sequence evolution, we need to consider two homologous DNA sequences which diverged from O, the common ancestor, t time units ago (Fig. 1). The number of substitutions per site between the two sequences is given by

K = 2t ~-~fi ~-~ri., = - 2 t ~'~firii , i

j4:i

(4)

i

where fi is the frequency of nucleotide i in the ancestral sequence; the factor 2 comes from the fact that K is the sum of the numbers of substitutions per site in the two lineages. One problem in estimating K is that K depends on the initial frequencies fi, which we do not know. However, if the nucleotide frequencies are stationary, that is, their expectations do not change with time, then the initial frequencies can be estimated by the average frequencies in the present sequences. Thus, the assumption of stationarity greatly simplifies the estimation problem. We consider the case where the substitution process is stationary and time reversible [7]; this is called the SR model. Time reversibility means that R is restricted by the following condition: q'l'iRij = 77"jR.i i

(5)

for any i :~ j, where 77"i is the equilibrium frequency of nucleotide i (i = l, 2, 3,4 for A, T, C and G). Therefore, the rate matrix R is given by

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167

162

A

T

C

G

A

rl I

"B'2v I

"7T3U2

7"/'4S 1

T

q'gl Vl

r22

"B'3$2

q'g4U3

C

q'/'l v2

"B'2s 2

r33

7"/-404

G

'7/'1SI

77"2U3

"B'3U4

r44

Since rii = - ~ j 4 : i rij, e.g. rl 1 = - ( 'B'2Ul + q'/'3u2+ "/7"4s1) and since 7ri + ~'2 + zr3 + 77"4= 1, the SR model is a nine-parameter model. The SR model includes the one-parameter model, the two-parameter model and several other models as special cases. Consider the two-lineage scheme (Fig. 1). Since time reversibility implies that the substitution process from the common ancestor O to sequences X and Y is equivalent to that from X through O to Y (or from Y through O to X), the transition probability matrix from X to Y is given by P(2t)

=

e 2tR

.

(6)

By spectral decomposition, the diagonal elements of R can be expressed as 4

rii = ~

UikVki/~k ,

(7)

k=l

where ,~k (k = 1,2, 3, 4) is the ith eigenvalue of R, one of which is zero, say "~-4 = 0; Uik is the ikth element of the eigenmatrix U and vki is the kith element of matrix V = U -j . By putting (7) into (4), we have 3

K = -2t~bkAk,

(8)

k=l

where the constants bk (k = 1,2, 3) are given by 4

bk = ~

7"TiUikVki.

(9)

i=1

In the following, we shall discuss how to derive the formula for estimating K under a specified substitution model. We will start with the one-parameter and two-parameter models, for which the eigenvalues of R can be obtained analytically. Then, we shall discuss the nine-parameter SR model, for which the eigenvalues of R can not be obtained analytically. Finally, we will show that all these methods can be extended to the case where the substitution rate varies among nucleotide sites.

W.-H. Li, X. G u / P h y s i c a A 221 (1995) 1 5 9 - 1 6 7

163

3. One-parameter method In the one-parameter model, the eigenvalues of R are given by ,ll = f12 = .t3 = - 4 u a n d /~4 = 0 and so bl = b2 = b3 = 1/4. From (8), the number of substitutions per site can be simplified as

K = 6ut.

(10)

From the eigenvalues and eigenmatrix of R, the transition probability is given by 1

P#(t) =

3 --4ut

4 + ~e 1 1 -4ut ~-~e

•

if i = j, i f / 4=j.

( 1I)

Now consider the two-lineage scheme (Fig. 1). Let l ( t ) be the probability that, at time t, the nucleotides at a given site in the two sequences are identical to each other. Since the probability that both sequences have nucleotide i is }--~,41 f k P ~ ( t ) , where fk is the frequency of nucleotide k at the ancestor O, I t is given by 4

l(t) = Z

fk { p 2 (t) + p22(t ) + P23(t ) + p24(t ) } .

(12)

k=l

From ( 1 1 ) and after some simplifications, one can show that, regardless of the initial frequencies fk, l ( t ) is given by

l ( t ) = at + a 3 e -8ut •

(13)

Note that the probability that two sequences are different at a site at time t is p = 1 - 1 (t). Thus,

p = ~ 1 - e-8"t).

(14)

Therefore K = 6ut can be estimated by K = - ~3 ln(1 - 4p) . ,

(15)

where p can be estimated by the observed proportion of different nucleotides between the two sequences.

4. Kimura's two-parameter method Under Kimura's two parameter model, the eigenvalues of R are A1 = A2 = --2(s + v), A3 = --4v and -~4 = 0 ; the constants are bl = b2 = 1/2 and b3 = 1/4. Thus, K is given by K = 2(s + 2v)t.

(16)

164

W.-H. Li, X. G u / P h y s i c a A 221 (1995) 1 5 9 - 1 6 7

The transition probability under this model is given by

Pq( t) =

¼( 1 -q- e -art -q- 2e-2(s+L')t),

for i=j,

¼( 1 +

for transitions,

e -4vt

--

2e-2(s+v)t),

¼( 1 - e-4Vt),

(17)

for transversions.

To estimate K, let p and q be the probabilities that, at a given site at time t, the two sequences differ by a transition and a transversion, respectively. One can show that

P =2~

fk(PkAPkG + PkTPa:),

(18)

k

and q = (1 - I ( t ) ) - p ,

where l ( t ) is defined by (12). So we have

1 --4(s+c)t 1 --8vt p = ~1 - ~e + ~e , 1 --8vt . q = 1 _ ~e

(19)

Thus, K = 2(s + 2 v ) t can be estimated by K = -½ ln(1 - 2p - q) - ¼1n(1 - 2q),

(20)

where p and q can be estimated from the two sequences.

5. The general stationary and time-reversible model By the above approach, explicit solutions for estimating K have been obtained up to six-parameter models [ 8 ]. However, under the general SR model, it is difficult to derive an analytical formula for K because the eigenvalues of R cannot be expressed in an analytical form. Fortunately, we can solve this problem as follows [7]. Let zk be the kth eigenvalue of P(2t). By matrix theory, (6) implies that Zk = e 2ta~ •

(21)

S i n c e ,~4 - - 0 and z4 = 1, there are only three non-trivial equations in (21). Eq. (6) also implies that R and P(2t) have the same eigenmatrix. Thus, the number of substitutions per site can be estimated by 3

K = - ~bklnzk,

(22)

k=l

where zk and bk can be estimated from sequence data (see below). It can be shown that, under the SR model, all eigenvalues Zk (or Ak) are real [7]. To estimate Zk and bk from sequence data, we must estimate the transition probability matrix P(2t) first. Let Jij be the expected frequency of sites where the nucleotide is i in sequence X and j in sequence Y. Let matrix J consist of Jij. It can be shown that matrix J is symmetric under the SR model. By Markovian properties, we have

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167

165

4

Ji) = Z q r k P k i ( t ) P k j ( t ) ,

i , j = 1 . . . . . 4.

(23)

k=l

By time reversibility, i.e. 7riPij(t) = wjPii(t), we have 4

4

Ji) = Z

TriPik( t ) Pkj( t) = 7ri Z

k=l

Pi~( t) Pkj( t ) = 7riPij( 2t) ,

(24)

k=l

where ~-]~4=1Pik(t)P~j(t) = Pij(2t) is a basic property of transition probabilities. Thus, (24) gives a simple method for estimating Pij(t) from sequence data Jij: (a) Count Nij, the number of sites at which the nucleotide is i in sequence X and is j in sequence Y, and then compute ]ij = N i j / N , where N is the number of nucleotides in the sequence. (b) Estimate the transition probability matrix P(2t) by

;,j

J'J

= 7=-, 77"i

i , j = l . . . . 4,

(25)

where ¢ri is the frequency of nucleotide i estimated by taking (simple) average between the two sequences. However, the estimated frequency Jo may not be equal to ]ji, i.e. the estimated matrix .] may not be symmetric. Li and Gu [9] proposed a method to test whether Jii - Jji is significantly different from zero. If not, the deviation from symmetry can be regarded as sampling effects, and the ijth element of j is replaced by l ( Jij -}- J j i ) .

Then, zk (k = 1 . . . . . 4) can be obtained by solving the characteristic equation d e t ( P zI) = 0, where P consists of/3/j and I is the identity matrix; the corresponding eigenmatrix U and its inverse matrix V are also obtained simultaneously by a standard algorithm. Thus, K can be estimated by (22) and the sampling variance can be approximately computed by the formula given by Gu and Li [7].

6. E s t i m a t i o n of K under variable rates

All the above methods for estimating K assume the same substitution rate for all sites. If this assumption is violated, K may be seriously underestimated, especially for divergent sequences. Because the substitution rate usually varies among sites in most genes, these methods need to be modified to take rate variation into account. There is evidence that the rate variation among sites follows a gamma distribution. This distribution is mathematically simple and is commonly used in the literature. Assume that the ijth element of the rate matrix R is expressed by rij = hiju, where the constant hi i represents the pattern of substitution rate and the random variable u varies among sites according to the following gamma distribution

&(u) = - - u ~ - l

F(a)

e -~u ,

(26)

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167

166

a/ft.

where the mean of u is given by ~ = First, consider the one-parameter model; in this case, = 1 for all i v~ j. From (14), the mean of p, i.e. the mean proportion of differences between two sequences, is given by

hij

~=fpqb(u)du=3{1-(l+~)-a}.

(27)

0

Therefore, the average number of substitutions per site K = 6~t can be estimated by

In the same manner, Jin and Nei [ 10] extended Kimura's two-parameter method [6] to the case where the rate varies according to a gamma distribution: K=~

°{2 ( 1 - 2 P - Q ) - ~ / " ÷ ( 1 - 2 Q )

-1/"-3

).

(29)

Finally, for the SR model, the average K is 3

K

=

,

_bk(zE

- 1),

(30)

k=l

where the eigenvalues zk and the constants bk can be estimated by the same approach as above for the SR model under the uniform rate model.

7. Concluding remarks Because the ancestral sequence is generally unknown, the most general model (i.e. 12 parameters) is difficult to apply so that some restrictions on the rate matrix R are necessary for developing useful methods for estimating K. As reviewed in this article, many methods have been developed for this purpose in the past two decades. As a general rule, methods based on more general models will have smaller estimation bias but larger sampling variances. Thus, when the sequence length (N) is large, more general methods are preferred. However, when N is small, simpler methods may be better. For example, Gu and Li's [7] simulation study suggested that, if N is less than 200 base pairs, the one-parameter method is on average better than the SR method; whereas if N is larger than 500, the SR method is on average better than the one-parameter method.

References 11 ] 121 13] 141

J. Felsenstein, Annu. Rev. Genet. 22 (1988) 521. C.I. Wu and W.H. Li, Proc. Natl. Acad. Sci. USA 82 (1985) 1741. A.L. Hughes and M. Nei, Proc. Natl. Acad. Sci. USA 86 (1989) 958. K.H. Wolfe, M. Gouy, Y. Yang, P.M. Sharp and W.H. Li, Proc. Natl. Acad. Sci. USA 86 (1989) 6201.

W.-H. Li, X. Gu/Physica A 221 (1995) 159-167

167

[5] T.H. Jukes and C.R. Cantor, in: Mammalian protein metabolism, H.N. Munro, ed. (Academic Press, New York, 1969) p. 21. 16] M. Kimura, J. Mol. Evol. 16 (1980) 111. 171 X. Gu and W.H. Li, Proc. Natl. Acad. Sci. USA (1995), submitted. 181 K. Tamura and M. Nei, Mol. Biol. Evol. 10 (1993) 512. [9] W.H. Li and X. Gu, in: Methods in Enzymology (1995), in press. I101 L. Jin and M. Nei, Mol. Biol. Evol. 7 (1990) 82.