Sunday, February 28, 2010

Jukes-Cantor (6)

I'm exploring the derivation of equations for the probability that a nucleotide either stays the same over a period of evolutionary time, or changes to be one of the other three. The Jukes-Cantor model is that all of these rates (X to Y) are the same, described by a parameter α that is the instantaneous rate.

One of my sources is Ziheng Yang's Computational Molecular Evolution. According to him (and wikipedia), this model can be represented by a substitution-rate matrix which is usually called Q:


     A     C     G     T
A -3λ λ λ λ
C λ -3λ λ λ
G λ λ -3λ λ
T λ λ λ -3λ


Yang uses λ in place of the often-used α. For this post we'll follow along with him.

As I posted about the other day, we need to have equations that will give us the probabilities for any time t. These should be recognizable:


PXX(t) = 1/4 + 3/4*e-4*λ*t
PXY(t) = 1/4 - 1/4*e-4*λ*t


Yang calls these p0(t) and p1(t).

The corresponding transition-probability matrix is designated P(t) and it is:


P(t):

p0(t) p1(t) p1(t) p1(t)
p1(t) p0(t) p1(t) p1(t)
p1(t) p1(t) p0(t) p1(t)
p1(t) p1(t) p1(t) p0(t)


And if we had P for one unit of evolutionary time, then we could do matrix multiplication (see here) to generate P for any time.

Now, Yang (and other sources) also says that:


P(t) = eQ(t)


And I've been trying to wrap my head around that, i.e. matrix exponentiation. Luckily, I've been getting into Gilbert Strang's MIT lectures on linear algebra (here). So I can appreciate what the inverse of a matrix is, such that:


Q Q-1 = I


where I is just the identity matrix of the same size and shape as Q. Now, suppose we can diagonalize Q. That is, we can find an invertible matrix U and a diagonal matrix Λ such that:


Q = U Λ U-1


A diagonal matrix means that the matrix has non-zero values only on the diagonal, e.g.:


λ1    0    0    0
0 λ2 0 0
0 0 λ3 0
0 0 0 λ4


And that means that Λ Λ is simply:


λ12    0    0    0
0 λ22 0 0
0 0 λ32 0
0 0 0 λ42


Look what happens when we do something like:


Q Q = U Λ U-1 U Λ U-1
= U Λ2 U-1


We only have to deal with Λ, and it's straightforward. Now for the real magic. I am reliably informed that this trick works for any function. For example:


eQ(t) = U e&Lambda(t) U-1


where we have:

e&Lambda(t) =

eλ1 0 0 0
0 eλ2 0 0
0 0 eλ3 0
0 0 0 eλ4


And, most impressive of all, U and U-1 are the right and left eigenvectors of Q, and the coefficients of Λ are its eigenvalues! We've looked at that in detail previously(here).

No comments: