Monday, September 9, 2013

Convert a tree with a mapped character to a tree with singleton nodes

The way that phytools stores a tree containing a mapped discrete character (for instance, created using the functions make.simmap or sim.history) is as a standard object of class "phylo", with an additional element ($maps) which is a list of the times spent in each state (in temporal sequence from the start to the end of the edge) along each edge.

Another way to store the same information would be as a modified "phylo" object in which each edge has one state, but so-called "singleton nodes" (nodes with only one descendant edge) are permitted. Converting from phytools standard format for mapped trees to this type of object was requested of me recently by Luke Mahler. The main reason for this was to be able to convert to an object of class "ouchtree". OUCH is a package that can be used to fit the multi-optimum Ornstein-Uhlenbeck model of Butler & King (2004). OUCH does not permit multiple regimes along a single branch in the tree (although this should be allowed theoretically by the model); however the OUCH function ape2ouch (which converts between an object of class "phylo" and an "ouchtree" object) accepts singleton nodes.

I just posted this function (map.to.singleton). Just so we can see that it's doing exactly what we thing it's doing, I have also added to phytools a new function plotTree.singletons which can do basic tree plotting of trees containing singleton nodes. Here's a quick demo of how the whole thing works:

> ## load phytools
> require(phytools)
> packageVersion("phytools")
[1] ‘0.3.46’
> ## first, let's get our tree with a mapped character
> tree<-pbtree(n=40,scale=1)
> Q<-matrix(c(-2,1,1,1,-2,1,1,1,-2),3,3)
> colnames(Q)<-rownames(Q)<-letters[1:3]
> tree<-sim.history(tree,Q,anc="a")
> tree

Phylogenetic tree with 40 tips and 39 internal nodes.

Tip labels:
        t22, t35, t36, t27, t28, t12, ...

Rooted; includes branch lengths.
> ## plot it
> plotSimmap(tree,lwd=2,pts=TRUE)
no colors provided. using the following legend:
      a        b        c
"black"    "red" "green3"
> ## ok, now let's convert to a tree with singletons & plot
> singles<-map.to.singleton(tree)
> singles

Phylogenetic tree with 40 tips and 72 internal nodes.

Tip labels:
        t22, t35, t36, t27, t28, t12, ...

Rooted; includes branch lengths.
> plotTree.singletons(singles)

That's pretty much it.

5 comments:

  1. Huh.

    Have you considered the opposite, a function that would convert a tree with singleton nodes and states for each edge, converted to a mapped character tree?

    I might have a use for such a thing... keeping track of ancestral taxa and ghost branches on a paleo-tree, for example.

    ReplyDelete
    Replies
    1. That shouldn't be too hard. If you have a specific example, then email me & I'll do it. Liam

      Delete
  2. Some reports from Luke suggest that this may depend on the latest non-CRAN ape (download here). I'm not going to roll back my version to check, but I will update if this is confirmed.

    ReplyDelete
  3. Just to follow up, that was indeed true, but the CRAN version of ape was updated the day after this post, so the CRAN version (3.0.10) works fine now.

    ReplyDelete
  4. Hello Liam,

    Thanks for adding this functionality. I trying to use the plotTree.singleton method in my scenario with the end goal of plotting ancestral state marginal probabilities as pie charts unto trees with singleton edges. To create these pie-charts, I'm using nodelabels and tiplabels functions (from ape). However, these functions don't seem to be supported on these singleton trees. Do you get errors if you try to call any of these two functions after calling plotTree.singletons in any of the above examples?

    In my scenario I'm running into unusual behaviour from the nodelabels function when I try to call these functions after plotting my tree with singletons edges. I can send a complete repro script with data and all, but a quick snippet of what I'm doing looks something like this:

    packageVersion("phytools")
    [1] ‘0.4.56’
    packageVersion("ape")
    [1] ‘3.3’

    mt = read.newick(text="...")
    plotTree.singletons(mt); #Looks good
    nodelabels(); #No errors, but the labels are not overlapping the ancestral nodes
    nodelabels(pie=marginals, piecol = cols); #ERROR
    Error in pie[i, ] : subscript out of bounds

    Thanks for taking a look at this. Do let me know if plotting pie charts at the ancestral nodes for singleton trees are possible.

    ReplyDelete