Monday, July 11, 2011

Writing a "simmap" tree to file

Ok, I just completed a very basic function to write a stochastic character mapped tree to file. The function is built on the general (if somewhat inelegant) algorithm for writing a "phylo" object to a Newick string that I gave in a previous post. This function is available from my R phylogenetics page (direct link to code here) and it will eventually join the "phytools" package.

Here's a quick example of potential usage. Let's start with our Newick tree from before:

> tree<-read.tree(text="(((2:1.55,((6:0.16,7:0.16):0.47,(4:0.22,5:0.22):0.41):0.92):0.24,1:1.79):0.17,(3:0.49,8:0.49):1.47);")

Now, let's simulate a discrete character on this tree using sim.char() from the "geiger" package:

> require(geiger)
> x<-sim.char(tree,model.matrix=list(matrix(c(-1,1,1,-1),2,2)), model="discrete")[,,1]
> x
2 6 7 4 5 1 3 8
1 2 2 1 2 2 1 1

Next, we can use make.simmap() to simulate a stochastic character map:

> mtree<-make.simmap(tree,x)

And we can plot it:

> cols<-c("blue","red"); names(cols)<-c(1,2)
> plotSimmap(mtree,cols,lwd=3,pts=F)

Finally, let's write the tree to a string:

> write.simmap(mtree)
[1] "(((2:{1,1.55},((6:{2,0.16},7:{2,0.16}):{1,0.01292597:2,0.13960021:1,0.22833548:2,0.08913835},(4:{1,0.22},5:{1,0.21640896:2,0.00359104}):{1,0.41}):{1,0.92}):{1,0.24},1:{1,0.0532546:2,0.59924269:1,0.1929242:2,0.94457851}):{1,0.17},(3:{1,0.24224632:2,0.09203948:1,0.15571419},8:{1,0.49}):{1,0.17465378:2,0.4633267:1,0.13014025:2,0.64147628:1,0.060403});"

Obviously, we can also write this to file. It is not built into write.simmap() yet, but we can just do:

> write(file="treefile.tre",write.simmap(mtree))

Or multiple trees:

> write(file="treefile.tre",write.simmap(make.simmap(tree,x)))
> write(file="treefile.tre",write.simmap(make.simmap(tree,x)), append=T)

Which we can then, of course, read in using read.simmap():

> trees<-read.simmap("treefile.tre")
Read 2 items
> plotSimmap(trees[[1]],cols)
> plotSimmap(trees[[2]],cols)


No comments:

Post a Comment