phytools has a function called `markChanges`

that can be
used to highly changes of different types on the tree. Various uses of
this function are documented
here;
however they focus on marking the changes from a single stochastic map
history. Alternatively, we might like to look at the *distribution*
of changes across a sample of stochastic maps.

This turns out to be really easy to do.

So, let's start with the following tree & data:

```
library(phytools)
tree
```

```
##
## Phylogenetic tree with 31 tips and 30 internal nodes.
##
## Tip labels:
## X1, A, B, C, D, E, ...
##
## Rooted; includes branch lengths.
```

```
x
```

```
## X1 A B C D E X2 X3 F G X4 H I J K L M N
## "a" "b" "b" "b" "b" "b" "a" "a" "b" "b" "b" "a" "a" "a" "a" "a" "a" "a"
## O P X5 Q R S T U V W X Y Z
## "a" "a" "a" "a" "b" "b" "a" "a" "a" "a" "a" "a" "a"
```

```
colors<-setNames(c("blue","red"),c("a","b"))
dotTree(tree,x,colors=colors)
```

We can do stochastic mapping as follows:

```
## first let's pick a model
fitER<-fitMk(tree,x,model="ER")
fitER
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b
## a -0.027553 0.027553
## b 0.027553 -0.027553
##
## Fitted (or set) value of pi:
## a b
## 0.5 0.5
##
## Log-likelihood: -10.956725
```

```
fitARD<-fitMk(tree,x,model="ER")
fitARD
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a b
## a -0.027553 0.027553
## b 0.027553 -0.027553
##
## Fitted (or set) value of pi:
## a b
## 0.5 0.5
##
## Log-likelihood: -10.956725
```

```
LR<--2*(logLik(fitER)-logLik(fitARD))
P<-pchisq(LR,df=1,lower.tail=FALSE)
P ## we cannot reject "ER"
```

```
## [1] 1
```

```
trees<-make.simmap(tree,x,nsim=100,model="ER")
```

```
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## a b
## a -0.02755304 0.02755304
## b 0.02755304 -0.02755304
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## a b
## 0.5 0.5
```

```
## Done.
```

```
trees
```

```
## 100 phylogenetic trees with mapped discrete characters
```

Now, let's make our plot. Unlike in the
example
referenced above, here I'm going to overlay the changes on on a
`dotTree`

plot, instead of on any single stochastic map. To make
the density of the marked changes a little bit easier to read, I'm going
to make them semi-transparent with the aide of the phytools internal
function `make.transparent`

(which I borrowed from something
Josef Uyeda shared with me a long
time ago):

```
dotTree(tree,x,colors=colors)
mt<-phytools:::make.transparent ## using the internal function
nulo<-sapply(trees,markChanges,sapply(colors,mt,0.4))
add.simmap.legend(colors=sapply(setNames(colors[2:1],c("a->b","b->a")),
mt,0.4),prompt=FALSE,x=0.5,y=30,vertical=TRUE)
```

If you wanted to go crazy, you could also add the posterior probabilities at nodes as follows:

```
obj<-summary(trees)
obj
```

```
## 100 trees with a mapped discrete character with states:
## a, b
##
## trees have 5.1 changes between states on average
##
## changes are of the following types:
## a,b b,a
## x->y 4.25 0.85
##
## mean total time spent in each state is:
## a b total
## raw 141.2496391 39.5391057 180.7887
## prop 0.7812966 0.2187034 1.0000
```

```
dotTree(tree,x,colors=colors)
nulo<-sapply(trees,markChanges,sapply(colors,mt,0.4))
add.simmap.legend(colors=sapply(setNames(colors[2:1],c("a->b","b->a")),
mt,0.4),prompt=FALSE,x=0.5,y=30,vertical=TRUE)
nodelabels(pie=obj$ace,cex=0.6,piecol=colors)
```

Although that could be a bit excessive.

In any case, it gives a good idea of where on the tree changes in state
most likely occured. For instance, it is very interesting to note that
changes `a->b`

most likely occurred independently on the
lineages leading to `X4`

and `(G,F)`

, rather than
once in their common ancestor followed by a reversion `b->a`

.

That's all.

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.