Tuesday, December 7, 2010

Models of sequence evolution 3

I wrote a couple of long posts about models of sequence evolution (here and here), in which (and it may not be obvious) we solved the problem of how to obtain:

P = eQt

for any t, by getting eigenvalues and eigenectors, using both R and Python. That's outstanding.

If you look carefully, you might notice something. There are just two terms in Q for the JC69 model:

-0.01 (on diagonal) and 0.003333 (off diagonal)

P has only two terms as well and they turn out to be (for t = 0.1):

0.999 and 0.000333

That is, the off-diagonal terms are apparently just multiplied by t, and naturally the diagonal is 1 - 3 times the other term. There was really no need for all the exponentiation of matrices, etc.

How does this happen? The formulas for the two terms in JC69 are:

PNN(t) = 1/4 + 3/4 exp(-4λt)
PNM(t) = 1/4 - 1/4 exp(-4λt)


But going back to the infinite series expansion for e (or just remembering the approximation for small x):

ex = 1 + x1/1! + x2/2! + ..
ex = 1 + x


So the exponential is just 1 - 4λt and we obtain:

PNN(t) = 1/4 + 3/4 (1 - 4λt)
= 1 - 3λt
PNM(t) = 1/4 - 1/4 (1 - 4λt)
= λt


Let's take another look at the formulas for these terms in the TN93 model. Look at only the first row of P. We first consider the two exponential terms:

e2 = exp(-βt)
e4 = exp[-(πYα1 + πRβ)t]

Approximating for small t:

e2 = 1 - βt
e4 = 1 - (πYα1 + πRβ)t

Look what happens to P13 and P14:


P13 = &piA(1 - e2)
= &piAβt

This is t times Q13. Exactly the same thing happens for P14. P12 is a bit more complicated, but it's just algebra:

P12 = πC + (πCπRY)e2 - (πCY)e4

The second term is:

πCπRY (1 - βt)
πCπRY - (πCπRY)βt

The third term is:

πCY (-1 + πYα1t + πRβt)

We can see that (considering both the second and the third term) there are two copies of: CπRY)βt with opposite sign, leaving:

P12 = πC + πCπRY - πCY + πCα1t

The last term is t times Q12. So I'm betting that the rest will go away. We have:

πC + πCπRY - πCY

but

πY + πR = 1
πR - 1 = -πY

so

πC + πCπRY - πCY
= πC + πCYR - 1)
= πC - πC
= 0



P12 = πCα1t = t Q12

How about that!

No comments: