Animating SIMMAP trees

It took me quite a while to get this working, but I finally did:

The tree topology is from Rüber et al. (2004).  This is what I used to generate the image:

animateSimmap = function(phy, interval = .02, numb = length(phy), name = "animation.gif", ...) {
  if(class(phy) != "multiPhylo") stop("object must be multiPhylo")
  saveGIF(for(i in 1:numb) {
            dev.hold()
            plotSimmap(phy[[i]], ...)
            Sys.sleep(interval)
          }, movie.name = name)
}

So this function takes in a multiPhylo object that has stochastic character mappings, the time interval between trees (I'm not sure if this does anything because saveGIF might just have a set interval already), the number of mappings to animate, the name of the file to create, and any arguments you want to give to the plotSimmap function.  In terms of the mappings shown above, they're meaningless because the branches aren't proportional to anything meaningful, but it shows how this function would work.

plotSimmap is from the package phytools, and saveGIF is from the package animation.  To use saveGIF, you also need either ImageMagick (which is what I used), GraphicsMagick, or LyX.  

SIMMAP Trees

Thanks to Liam Revell, we can now produce simulated discrete character mappings within R.  The function make.simmap uses ace, from the package ape, to fit the model to use for the simulations.  Thus, the proportion at which a character state appears at a node among many iterations of simulations should be roughly equivalent to the likelihood of that state as estimated in ace.  I have heard this from many people, but I wanted to be able to summarize the actual simulated states.  In order to do this, I took a look at the components of a SIMMAP tree.  I'll illustrate this with a mock example:

We can take a look at what components this object has:

> names(exampleSimmap)
[1] "edge"        "edge.length" "tip.label"   "Nnode"       "maps"       
[6] "mapped.edge"

The first components are the same as any phylo object.  So maps and mapped.edge are what make a SIMMAP tree special.  Let's take a look (the middle elements removed to save space):

> exampleSimmap$maps
[[1]]
        1         0 
0.0958356 0.3197055 
[[2]]
        0         1 
0.2371619 0.1659281 
[[3]]
         1 
0.01380231 
...
[[17]]
        1 
0.1164989 
[[18]]
        1         0 
0.1161601 0.1769559 

This is precisely what we need!  Each element of exampleSimmap$maps represents a single branch, and the values represents the length of time that branch spends in each state, in this case 0 or 1.  That means we can simply take each branch's starting value (whatever is the name of the first element of that branch), and that is the value at the node where the branch starts.  Let's see if we can find this.

> exampleSimmap$edge
      [,1] [,2]
 [1,]   11   12
 [2,]   12   19
 [3,]   19    1
 [4,]   19    2
 [5,]   12   13
 [6,]   13   17
 [7,]   17    3
 [8,]   17    4
 [9,]   13   18
[10,]   18    5
[11,]   18    6
[12,]   11   14
[13,]   14   15
[14,]   15    7
[15,]   15   16
[16,]   16    8
[17,]   16    9
[18,]   14   10

The element named edge gives us the starting and ending node for each of the 18 edges in our tree.  That means we can use to figure out which node corresponds with which state.  Here I've written a function that takes in a SIMMAP tree and returns a named vector where the values are the node states and the names are the nodes.

mappedNode = function(phy) {
  # phy must be a SIMMAP tree
  nodes = phy$edge[,1] # this gives us the starting node for all edges
  map = sapply(phy$maps, function(x) attr(x, "names")[1]) # this gives us the starting value of each branch
  df = unique(data.frame(nodes = nodes, map = map)) # here we're removing the repeated values as interior nodes will have multiple branches
  mapping = df$map
  names(mapping) = df$nodes # naming the vector with node names
  mapping
}

Now we can run this function over all of the simulated mappings:

mappings = sapply(simmapTrees, function(x) mappedNode(x))
# change the 0 and 1 to numeric
nummaps = as.data.frame(sapply(1:length(simmapTrees), function(x) as.numeric(mappings[,x])))
# make sure the row names correspond to node names
rownames(nummaps) = rownames(mappings)
# get the number of simulated trees with  node state of 1
sums = sapply(rownames(nummaps), function(x) sum(nummaps[x,]))
# change that to a frequency
freq = sums/length(simmapTrees)
# plot
plot(exampleTree, label.offset = .05)
nodelabels(pie = freq, cex = .65, node = as.numeric(names(freq)))
tiplabels(pie = discTrait, cex = .65)
Looks fairly reasonable.  Now let's estimate the ancestral character using ace:

MLACE = ace(discTrait, exampleTree, type = "discrete", model = "SYM")

plot(exampleTree, label.offset = .05)
nodelabels(pie = 1 - MLACE$lik.anc, cex = .65)
tiplabels(pie = discTrait, cex = .65)

Pretty close, but not identical.  Maybe if we do more simulated mappings (I did 1000 here), they will start to look more similar?

Regular Expressions

My first experience with regular expressions came from Python for Dummies.  It looked particularly relevant to the specific task I was working on, scraping specific bits of information from fishbase.  When my advisor, Peter Wainwright, first approached me about this, I didn't know where to begin, so I went to two people with experience in these sorts of tasks:  Bob Thomson and Carl Boettiger.

Bob's suggestions were to download each individual fish's HTML file using a short bash script, then use something like Python or Perl to extract relevant bits.  With over 30,000 species, just downloading the HTML files took quite a long time.  But with no background (at the time) in Python or Perl, I turned to Carl, who suggested using R.  He quickly wrote a package, rfishbase that allows you to access information from the XML files on fishbase through R.  Although the XML files don't have all of the information available on the HTML files, they still have quite a lot.

My reason for this post, though, is because of a task my lab mate, Patrick, wished to accomplish using the data he accessed using rfishbase.  Looking at a character vector containing information of interest, he wanted to get all of the reference numbers within that vector.  An example of an element might be something like this:

"Occurs mainly over rocky and muddy bottoms.  Uncommon around coral reefs.  Usually rests on the bottom (Ref. 9710).  Juveniles may be found in shallow water, but adults are usually taken from depths of 70-330 m (Ref. 13442).  Reptant and natant decapods were the main food items throughout the year (Ref 59311).  Feeds on a wide variety of fishes and invertebrates."

Given this, he would want the numbers 9710, 13442, and 59311.  Even in this one example, you can see that they are not always consistent:  the first two have a period while the third doesn't.  And there are even things like this:

"Common species.  Free-living.  Assumed to feed on small invertebrates and fish (Ref.  4741, 34024).  Feed on small bottom animals (Ref. 35388)."

Notice the many spaces before the first reference and having two numbers.  Or this:

"Occurs in various inshore habitats (Ref. 9800). Feeds on benthic invertebrates and fish (Ref. 11889). Also Ref. 43081."

This one doesn't even have parentheses around the last one.  So the first thing I did was to find every instance of "Ref" followed by an optional period, any number of spaces, and a run of any number of numbers, commas, and spaces.

ref = regmatches(matches, gregexpr("Ref\\.? *[0-9 ,]*", matches))

where matches is the character vector with all of the information we're looking at.  This is modified from here.  Next, I removed all characters other than digits or commas and used strsplit to separate individual reference numbers.

refs = sapply(ref, function(x) unlist(strsplit(gsub("[^0-9,]", "", x), ",")))

You end up with a list the same length as the original character vector, and every element is a character vector of all of the reference numbers.  From here, you can go in and find all of the unique values to find all of the references you need.

Geiger Bug

Big thanks to Luke Harmon for supplying me with updated code that fixes a small bug in the geiger package!  For those of you who have had issues with fitDiscrete in the past, you can email Luke for a fix.

The issue that arises is when you use either the symmetrical model (model = "SYM") or the all rates different model (model = "ARD").  You will get an error message looking like this after waiting for some time for the likelihood optimization to finish:

Finding the maximum likelihood solution
[0        50      100]
[....................]
Error in getQ(exp(out$par), nRateCats, model) :
  You must supply the correct number of rate categories.

But you don't actually control what goes into getQ through the arguments you give to fitDiscrete.  So if you encounter this issue using either of the above models, be sure to email Luke to get the updated code.  A revamped version of geiger is on the way, so hopefully this issue won't be around for much longer!

The one thing Luke cautioned me about using these two models is that they can both quickly become parameter-rich.  The number of rate parameters for the symmetric model is n * (n-1) / 2, while for the all rates different model, the number of rate parameters is n * (n-1) where n is the number of discrete states.  For example, if I have five discrete states, I would have ten parameters for the symmetric model and twenty parameters for the all rates different model.  I need a lot of data to be estimating so many parameters!

I did try out his updated code, and it works perfectly fine.  It still does take a while (I'm sure the speed depends on the size of your data set, the shape of your phylogeny, your computer's specs...), but I hear that there is already a faster version of the code if you ask Luke for it.  It would definitely be interesting to test out how much of a difference there is and what they changed to make it faster.

If all you are trying to estimate is a single rate model (model = "ER", which stands for equal rates), then there is no need to use an updated version of the function.  The old version will work just fine.  The problem with the original was that the calculation of the rate categories occurred twice:  once to get the number of rate categories, and another within getQ to test whether the number of rate categories was correct.  So for example, with five discrete states nRateCats in the above will equal ten for the symmetric model.  But you will get an error to supply the correct number of rate categories because within getQ, it compares the number of parameters to nRateCats * (nRateCats-1) / 2, which would compare ten to forty-five.

Creating Trees

I told Luke Mahler about how I'm interested in practicing analyses in R with my particular group of interest, even though I haven't collected data on more than a few species and I don't have a tree.  He told me creating a tree from previously published studies is easy:  just import the tree into R in newick format! In the geiger package in R, there is a function read.tree that does just that.

The tricky part about doing this is creating the newick-format tree itself.  It involves a lot of parentheses, colons, and parentheses:

exampleTree = read.tree(text = "(((A:1, B:1):1, C:2):1, D:3);")

will give you

Since I was typing out the newick-format tree directly within read.tree, I needed to make sure it went into the text argument instead of the file argument.  As for the newick format itself, the tip names are followed by the immediately subtending branch length:  hence A:1 or D:3.  The nodes also need to be provided with branch lengths, which is the :1 following each clade designated in parentheses. Finally, you can't forget the semi-colon at the end.  It is very easy to get lost in a sea of parentheses, especially for large trees, so inevitably I needed to go back to make changes.  Testing small clades at a time makes this a little easier.

In my particular tree, I had polytomies, which I designated with 0-length internal branches.  So if in the above tree, I actually don't have any information about the interrelationships of A, B, and C, I can use this:

polytomyTree = read.tree(text = "(((A:1, B:1):0, C:1):1, D:2);")

to get

RStudio

I was first introduced to RStudio when I took STA141 with Duncan Temple Lang.  At the time, I was using a PC, and it made things much simpler.  A few of the features that I love about it:

1) Color scheme for scripts.  I can see at a glance what's a comment, and the parentheses pop out at you.  This was something I envied about Macs when I used to use my PC exclusively.  I do wish you could personalize it a little more, but the schemes they have are great.  I personally like Cobalt.

2) Balancing parentheses and quotations.  Yes, this can get obnoxious if you are tweaking existing code rather than starting from scratch, since if you try to type an end quote, it will be interpreted as another beginning quote.  But I find being able to see what's matching your current end parenthesis incredibly helpful.  Besides, if you don't like it, you can easily turn it off in the preferences.

3) Everything in one window.  This wasn't a step up from my PC version since the normal R console keeps everything in one window, but it's so nice to have everything (script, console, history, help pages, graphics) in one window.

Color-coding

I've been helping some fellow students here at the workshop with R.  One of the things I did was to create a vector of colors to label the tiplabels.  Basically, the structure of the data was a data frame with the ID that roughly corresponds with the tip labels in the first column and the actual species name in the second column.  Here being a toy example:

> species
     ID species
1  A123   speca
2 SA123   speca
3  c123   specb

However, the tip labels of the tree don't exactly correspond to the IDs:

> tree$tip.label
[1] "SA123" "a123"  "c123"

Time for some regular expressions!

> species[,1] = gsub('^A(.*)', 'a\\1', species[,1])
> species[,1]
[1] "a123"  "SA123" "c123"

I'll break down the gsub function call a little bit.  The first argument to gsub is the pattern.  The ^ means the beginning, then A.  () names a group, and since it is the first set of (), this group is 1.  Within 1, there is a .*.  The . indicates any character, and the * indicates any number of those characters.  So essentially, this pattern is looking for a capital A at the beginning of a string followed by anything or nothing.

The second argument is the replacement.  The sub in gsub stands for substitution.  So it will take the entire pattern and replace it with whatever the replacement is.  If nothing matches, then nothing happens.  The \\1 refers to the named group, 1.  The \\ escapes, so that gsub will replace with whatever it finds in group 1, not with the number 1.  So essentially, it will match any string that starts with a capital A followed by a group 1 that can be anything, and replace it with a lowercase a followed by the characters in group 1.  The third argument is simply the vector you want to perform this on.

Now that I've gotten the ID to match the tip labels, I'm going to reorder the data frame to match the order of the tip labels.

> rownames(species) = species$ID

> species = species[tree$tip.label,]
> species
         ID species
SA123 SA123   speca
a123   a123   speca
c123   c123   specb

Now, all I have to do is change species$species into a factor and convert it to a numeric!