Transition probability of the HKY model

2005-10-04 shimo

Some identities and transformations which will be used later

Off[General::spell]

id1 =  {(Subscript[π, A]) + (Subscript[π, G]) + (Subscript[π, C]) + (Subscript[π, T]) == 
            1, (Subscript[π, Y]) + (Subscript[π, R]) == 1};

tr1 = {(Subscript[π, C]) + (Subscript[π, T]) -> (Subscript[π, Y]), (Subscript[π, A]) + (Subscript[π, G]) -> (Subscript[π, R])};

The instantaneous transition matrix

hkyA=(({{-β (Subscript[π, C]) - α (Subscript[π, G]) - β (Subscript[π, T]), α (Subscript[π, A]), β (Subscript[π, A]), β (Subscript[π, A])}, {α (Subscript[π, G]), -α (Subscript[π, A]) - β (Subscript[π, C]) - β (Subscript[π, T]), β (Subscript[π, G]), β (Subscript[π, G])}, {β (Subscript[π, C]), β (Subscript[π, C]), -β (Subscript[π, A]) - β (Subscript[π, G]) - α (Subscript[π, T]), α (Subscript[π, C])}, {β (Subscript[π, T]), β (Subscript[π, T]), α (Subscript[π, T]), -β (Subscript[π, A]) - α (Subscript[π, C]) - β (Subscript[π, G])}}));

which satisfies

{1, 1, 1, 1} . hkyA

The expected rate of change

The expected rate of transitions are obtained as

hkyE = 
      hkyA . DiagonalMatrix[{(Subscript[π, A]), (Subscript[π, G]), (Subscript[π, C]), (Subscript[π, T])}];  
    hkyE // MatrixForm

Simplify[{1, 1, 1, 1} . hkyE, id1]

The expected rate of total change is

hkyT = Expand[-Sum[hkyE[[i, i]], {i, 1, 4}]]

which may be simplifed as

hkyT = Collect[hkyT, {α, β}, Simplify] /. tr1

Let r = α/β. Then, α and β are defined by giving a constaint that hkyT=1.

Solve[(hkyT /. {α -> r β}) == 
        1, β]

The eigenvalues

foo = Eigensystem[hkyA]

The eigenvalues

foo[[1]] // TableForm

val = FullSimplify[foo[[1]], id1]

The eigenvectors

Transpose[foo[[2]]] // MatrixForm

hkyV = 
        Simplify[
          Transpose[foo[[2]]] . 
            DiagonalMatrix[{(Subscript[π, T]), 1, 
                1, ((Subscript[π, A]) + (Subscript[π, G])) (Subscript[π, T])}], id1];

hkyV // MatrixForm

This matrix consists of columns of eigenvectors

Check the eigenvectors

Define a diagonal matrix of eigenvalues

hkyL  =  DiagonalMatrix[val];

hkyL // MatrixForm

check if A*V = V*L

Simplify[hkyA . hkyV - hkyV . hkyL, id1]

Inverse matrix of hkyV

hkyIV = Simplify[Inverse[hkyV], id1];

hkyIV // MatrixForm

check if V*IV = Identity matrix

Simplify[hkyV . hkyIV, id1]

check if  V*L*IV=A

Simplify[hkyV . hkyL . hkyIV  -  hkyA, id1]

Transition probability matrix P(t)=Exp[A t]

Define Exp[L t]

hkyExpLt  =  DiagonalMatrix[Exp[val t]];

hkyExpLt // MatrixForm

Define Exp[A t] =  Exp[V*(L t)*IV] = V*Exp[L t]*IV

hkyExpAt = 
          Simplify[hkyV . hkyExpLt . hkyIV , id1];;

hkyExpAt // MatrixForm

Look at the each term of the transition probability matrix

Prob(C|A,t)

Transition probability from A to C (transversion) is

Pcat = hkyExpAt[[3, 1]]

Prob(G|A,t)

Transition probability from A to G (transition) is

Pgat = hkyExpAt[[2, 1]]

Try simplification of the above

Pgat = 
      Collect[Expand[hkyExpAt[[2, 1]] /. tr1 /. (Subscript[π, Y]) -> 
              1 - (Subscript[π, R])], {(Subscript[π, G]), (E^(-t β))}, Simplify[#] &]

The transition probability Prob(G|A,t) given in eq.(13.10) at page 203 of Felsenstein (2003) is

joePgat =  
        Exp[-β t] (1 - 
                Exp[-(Subscript[α, R]) t]) ((Subscript[π, G])/(Subscript[π, R])) + (1 - 
                Exp[-β t]) (Subscript[π, G]);

in which (Subscript[α, R]) for HKY model is defined as

hkytr1 = {(Subscript[α, R]) ->  (α - β) (Subscript[π, R]), (Subscript[α, Y]) ->  (α - β) (Subscript[π, Y])};

So, the Pgat for HKY model should be

joePgat /. hkytr1

Check if Pgat = joePgat

Simplify[Pgat  -  (joePgat /. hkytr1)]

Prob(T|C,t)

Transition probability from C to T (transition) is

Ptct = hkyExpAt[[4, 3]]

Try simplification of the above

Ptct = 
      Collect[Expand[hkyExpAt[[4, 3]] /. 
              tr1 /. (Subscript[π, C]) -> (Subscript[π, Y]) - (Subscript[π, T])], {(Subscript[π, T]), (E^(-t β))}]

\(Check\ if\ \ Ptct\ is\ obtained\ from\ Pgat\ by\ replacing\ π\_G -> π\_T, π\_R -> π\_Y\)

Simplify[
      Ptct  -  (Pgat /. {(Subscript[π, G]) -> (Subscript[π, T]), (Subscript[π, R]) -> (Subscript[π, Y])})]