BM simulation with a trend

Dave Bapst just posted a query to the R-sig-phylo about simulating Brownian evolution with a trend. I take this to mean that evolutionary changes along branches in the tree are drawn from a random normal distribution with variances sig2*tree$edge.length and means mu*tree$edge.length. [Initially, in my first response to the query, I did not multiply mu by tree$edge.length; but that doesn't make any sense.] I have now added this capability to v0.2 of my function fastBM() which can be downloaded from my R phylogenetics page here.

After loading the function from source, we can just type:

> x<-fastBM(tree,a=0,mu=0.5,sig2=2,internal=FALSE)

at the command prompt. This will return a n x 1 vector of character values (in x for the n species in tree.

To test out this function, let's try simulating data and then fitting a "trend" model using {geiger}'s fitContinuous() function. [Obviously, since these are stochastic simulations, individual results will vary.]

> tree<-rtree(500) # first simulate stochastic tree
> x<-fastBM(tree,mu=0.5,sig2=2) # now data on the tree
> fitContinuous(phy=tree,data=x,model="trend") # now fit model
Fitting trend model:
$Trait1
$Trait1$lnl
[1] -900.9902
$Trait1$mean
[1] 1.087359
$Trait1$beta
[1] 1.87744
$Trait1$mu
[1] 0.4968745
$Trait1$aic
[1] 1807.980
$Trait1$aicc
[1] 1808.029
$Trait1$k
[1] 3

Seems to work!

Note that Gene Hunt also pointed out (correctly, of course) that this can be done using the {MASS} function mvrnorm(). To see the discussion thread on this topic, check out the archived R-sig-phylo emails here.

New read.simmap() function

OK, I have an "alpha" version of a new read.simmap() function that can read in a stochastic character mapped (SIMMAP v1.0) style tree from file in "phylip" format and retains the ordering of the character states along the edges of the tree. Although it is incomplete (by definition), I have posted this function on my R phylogenetics development page here (it is v0.2.alpha of read.simmap()).

The way I decided to store the information about the order and time spent in each state on every edge was in a list of named vectors called $maps. For backward compatibility with my other functions using SIMMAP format trees (e.g., brownie.lite()) I also compute the matrix $mapped.edge containing the summed times spent in each state along each edge of the tree.

So, for the following phylogeny:

((A:{aqua,0.3:terr,0.4:aqua,0.3},B:{aqua,1.0}):{aqua,1.0},(C:{aqua,0.25:terr,0.75},D:{aqua,1.0}):{aqua,1.0});

the modified "phylo" object returned by my function has the following structure:

> temp<-"((A:{aqua,0.3:terr,0.4:aqua,0.3},B:{aqua,1.0}):{aqua,1.0},(C:{aqua,0.25:terr,0.75},D:{aqua,1.0}):{aqua,1.0});"
> tree<-read.simmap.alpha(text=temp)
> tree$Nnode
[1] 3
> tree$tip.label
[1] "A" "B" "C" "D"
> tree$edge
[,1] [,2]
[1,] 5 6
[2,] 6 1
[3,] 6 2
[4,] 5 7
[5,] 7 3
[6,] 7 4
> tree$edge.length
[1] 1 1 1 1 1 1
> tree$maps
[[1]]
aqua
1
[[2]]
aqua terr aqua
0.3 0.4 0.3
[[3]]
aqua
1
[[4]]
aqua
1
[[5]]
aqua terr
0.25 0.75
[[6]]
aqua
1
> tree$mapped.edge
state
edge aqua terr
5,6 1.00 0.00
6,1 0.60 0.40
6,2 1.00 0.00
5,7 1.00 0.00
7,3 0.25 0.75
7,4 1.00 0.00

Next, I will add the capacity to read "nexus" format SIMMAP v1.0 style trees.

Adding branch lengths

In a few different prior posts (1,2,3) I described the development of a simple R function, read.newick(), that is available from my development page here. This function just reads a single Newick style tree from a character string or from file, and then creates an {ape} "phylo" object in memory. The purpose of developing this function is to add the capability of reading stochastic character map information from within the Newick tree (a problem that I tackled a very different way in a prior post). I've also provided the code online so that readers can see the R implementation of our tree reading algorithm presented in a prior blog post.

Today I just added the reading of branch lengths to my previous function. This seems to work fine, and it is easy to see how I did this by reviewing my code - to which I have also added some comments.

Branch lengths in a Newick style tree are given after a full colon following either a label or a ")" right bracket. Thus, the branch lengths for our simplified primate tree might be: (((Human:7,Chimp:7):3,Gorilla:10):15,Monkey:25);, with branch lengths in millions of years (more or less). When our parser encounters the character ":", it just needs to assign the real number following the ":" to the edge preceding the node we are on. Otherwise, we just follow our algorithm as previously described.

The only amendment to the structure of the tree in memory is a vector, $edge.length, containing the length of each branch where the order corresponds to the order of the edges in $edge. Thus, for the tree (((Human:7,Chimp:7):3,Gorilla:10):15,Monkey:25);, we would get the following:

> temp<-"(((Human:7,Chimp:7):3,Gorilla:10):15,Monkey:25);"
> phy<-read.newick(text=temp) # or use read.tree()
> phy$Nnode
[1] 3
> phy$tip.label
[1] "Human" "Chimp" "Gorilla" "Monkey"
> phy$edge
[,1] [,2]
[1,] 5 6
[2,] 6 7
[3,] 7 1
[4,] 7 2
[5,] 6 3
[6,] 5 4
> phy$edge.length
[1] 15 3 7 7 10 25

Phylomorphospace updated

Luke Mahler identified an error in the node coloring option of my phylomorphospace() function. Basically, the function was assigning the colors of the daughter node to each parent node every time that a new edge was plotted. This has now been fixed and the updated version is available from my R-phylogenetics beta version page. Thanks Luke.

Reading a phylogeny from Newick (cont.)

In the last couple of posts, I have discussed a couple of basics regarding the input and storage in memory of phylogenetic trees in R. In particular, in the first post, I described the structure of the {ape} "phylo" object; and in the second, I described a standard algorithm (not of my own invention, but originally conveyed to me by Luke Harmon) for parsing a string in Newick format into a tree consisting of nodes and edges. The ultimate goal is not to read Newick trees into R (which read.tree() does quite well already), but to use this function as the basis for reading in other information from file - such as (for instance) the stochastic character map of a binary or multistate character.

Now our task is to translate the algorithm described in Saturday's post into computer code which will create the structure described in Friday's post. I have done this, and today I posted this function to my R-phylogenetics development page here; however, I remind the reader that this function is really provided for development/experimentation purposes only as it does nothing that is not already implemented in the {ape} function read.tree().

Ok, the first thing we need to do is figure out how many internal nodes and tips we have in the tree. We'd like to do this a priori, because that way we have two important pieces of information before starting: 1) how many rows to put in our matrix $edge; and (more importantly) 2) what number to assign to the root. The way I have done this is by taking one pass through the Newick string and counting the number of commas and the number of right brackets. It is my belief [readers, please correct me if this is wrong!] that the number of commas + 1 will always correspond to the number of tips in the tree; and that the number of right brackets will always correspond to the number of internal nodes.

Once we have figured out how many tips and internal nodes there is, it is straightforward to count the number of edges in the tree. This will just be the sum of the number of internal nodes + the number of tips - 1. The logic for this is as follows: every tip and internal node have an edge leading to them, except for the root node in the tree (hence the -1).

Now with these preliminaries completed, we can proceed to our algorithm for parsing the Newick string. Here, I start with a "while(tree[i]!=";")" loop that will terminate when the end character of our Newick style tree (";") is reached. Inside this loop, we can just program our character-by-character parser to build the tree based on the algorithm of my previous post and animation. I programmed this function more or less as illustrated, with a couple of caveats: 1) I never actually visit terminal nodes (I just create them as daughters, but then never update the status of "current node" to the daughter); and 2) in progressing backwards through the tree in response to character ")", the structure of the "phylo" object prevents easy backward movement through the tree. For this, I used the {base} function match() to match my current node's parent to the row number of prior nodes in the matrix edge.

Once all the elements of the tree are created, I just join them together in a named list using the function list(); and then I assign this list the class "phylo". Finally, I return the list to the R environment - and we're done!

Building a tree in memory

Ok, as discussed in my last post, our basic goal is to read a Newick style tree [for example, the tree (((Human,Chimp),Gorilla),Monkey);] into computer memory. We would like to do this by parsing the character string, element-by-element, to build a representation of the phylogeny that consists of nodes and edge. Of course this functionality already exists in the form of the {ape} function read.tree(); however, if we can develop a parsing function, we can then adapt this function to simultaneously read in other kinds of information - for instance, the order and time spent in each character in a phylogenetic stochastic character map.

Before showing how this can be done to create the {ape} "phylo" object discussed in the last post, I thought it would be a good idea to present the basic algorithm that is used to parse a Newick tree. I'm not sure if this handy set of rules has ever been published, but it was conveyed to me verbally by Luke Harmon a number of years ago.

The basic algorithm reads the Newick tree "character-by-character" (actually, we will read node labels as if a single character), with the following rule.

Starting with the root node, if we encounter:
1) "(" -> we will create left daughter node;
2) "," -> we will go back to the parent node and create right daughter;
3) ")" -> we will go back to the parent;
4) ";" -> we will finish; and, finally
5) if we encounter a character string [not including special characters listed above], then we label the node with that string.

To be a little clearer, I made a simple animation for the (((Human,Chimp),Gorilla),Monkey); tree from last time. The animation goes through each character of the Newick string at the top of the screen, and creates nodes and edges according the rules above. Enjoy.

Reading trees

In a comment on one of my earlier posts Carl Boettiger remarked that my read.simmap() R function (to read stochastic character mapped phylogenies) destroys the ordering of states along each branch and retains only the time spent in each state on each edge. This ordering is unimportant for analyses based around an undirected Brownian motion model - such as, for instance, the model of O'Meara et al. (2006); however, ordering could be important for other models, such as the Ornstein-Uhlenbeck model (e.g., Butler & King 2004). Carl is a sharp guy and his point is well taken, so now that I have finally settled into my new digs here (not-so-subtle advertisement - grad students welcome!), I decided to try to tackle this problem today.

The first thing I realized was that it would not be practical to use the same "work-around" algorithm that I used for v0.1 of read.simmap() and described here. I was going to have to figure out how to read phylogenies into memory in R without being able to take advantage of the {ape} function read.tree()!

To convey what this problem involves, it is first important to understand how a phylogenetic tree is stored in memory in R by {ape}. read.tree() stores the phylogeny as list of class "phylo" with 3 basic elements, as follows:

$Nnode : this is an integer containing the number of internal nodes in the tree;
$tip.label : this is a vector containing the labels for all the tips in the tree; and finally
$edge : this is a matrix of dimensions ($Nnode+length(tip.label)-1 x 2) which contains the node numbers for all the internal parent (left column) and daughter (right column) nodes in the tree.

By convention in {ape}, terminal nodes are numbered 1:n for n species, where the numbers correspond to the order of the tip labels in $tip.label; and internal nodes are labeled (n+1):($Nnode+length(tip.label)) starting with the root node of the tree.

[Note that this is really different from how phylogeny storage is typically programmed in C; where we would use a network of linked custom data types called "structures. An example of this can be seen on my website, here.]

To give you an idea of what this would look like then, the following tree (see figure) in Newick format:

>tree<-read.tree(text="(((Human,Chimp),Gorilla),Monkey);")

would have the following elements:

> tree$Nnode
[1] 3
> tree$tip.label
[1] "Human" "Chimp" "Gorilla" "Monkey"
> tree$edge
[,1] [,2]
[1,] 5 6
[2,] 6 7
[3,] 7 1
[4,] 7 2
[5,] 6 3
[6,] 5 4


More on this very soon. . . . .