Re: Two more questions

[Nick Matzke]

On Tue, Sep 12, 2017 at 2:01 AM, 张睿 <xiaor...@gmail.com> wrote:

Sorry, I have got two more questions:

How could I rotate a branch of phylogenetic tree?

The APE function for this is "rotate", I think.

 

And is there a function similar with “oder nodes decreasing / increasing” in FigTree?

I think the APE function is ladderize().

 

Thanks!

However, you cannot rotate nodes, stick that tree into an already-done BioGeoBEARS analysis, and then expect the plotting to come out OK.  This is because the state/range probabilities are stored in a specific order, and the order will change if you change the tree.  

It wouldn't be mind-bogglingly difficult to write some code to re-order the stored state/range probabilities to match some reordered tree -- you could do it with the BioGeoBEARS prt() function, which prints the tree to a table.  Make sure to set get_tipnames=TRUE to get the alphabetical list of the species names in each clade on the tree, this is a unique list that can be matched between the new tree and old tree with the match() function.

Easiest option: alter your tree in FigTree so it looks like you want, save that to Newick, use the new Newick file for re-running the BioGeoBEARS analysis with the new tree.

Note: I think FigTree displays trees "flipped" vertically, compared to how the identical Newick file will display in R.

Harder option: alter your tree in FigTree or R so it looks like you want, save that to Newick, use prt() and match() to re-order all of the probabilities tables in the saved BioGeoBEARS objects to match the new tree, give the BioGeoBEARS objects the new Newick file, then plot the old BioGeoBEARS results on the new tree.

Cheers,

Nick

PS: The links below will help -- but it will take some serious study of them:

Table of probabilities of states/ranges at each node

http://phylo.wikidot.com/example-biogeobears-scripts#node_probabilities

Fwd: state/range probabilities at each node in BioGeoBEARS...and other google group posts...

https://groups.google.com/forum/#!topic/biogeobears/eGzLRPg10y0

And see especially:

Some code for working with phylo objects in R http://luisdva.github.io/edge-lengths/  cc @jeff_o_hanson @NickJMatzke

https://twitter.com/LuisDVerde/status/836361904645943297

And this post on the use of prt():

(pasting the whole thing below for future Google-ability)

Making sense of trees

http://luisdva.github.io/rstats/edge-lengths/

==========================================

Making sense of trees

In the words of Nick Matzke:

R’s tree structure is pretty non-intuitive, compared to the “a tree is a collection of node objects” structure that is taught in most phylogenetics courses and used in e.g. C++ software. It’s done this way because R likes lists, and doesn’t like objects.

The structure of phylogenetic trees in R (namely APE’s phylo objects) can be confusing because much of the information needed to describe how different nodes and branches make up a tree is implicit.

I learned this the hard way during the early days of my PhD research. My thesis work revolved around phylogenetic lineage ages, so one of my very first tasks after finding a suitable tree was to get a table with species and their ‘edge lengths’ as a way to (more or less) represent their evolutionary ages. At this time I was very new to R and I lost months trying to figure this out on my own. After going over the APE book I slowly realized that lineage age information is contained within the phylo objects and that it can be extracted, but only if you ask nicely.

Phylo objects have a vector of edge lengths for every node, but I was unable to figure out how to get a list of tips and internal nodes that would match up with the vector of edge lengths. After much struggle, I figured out that the branching.times() function computes the distance from each node to the tips, producing a named vector that corresponds to the nodes in the tree. I also learned that the mrca() function produces a matrix of node numbers corresponding to the most recent common ancestors between pairs of tips or nodes. The next challenge would be to put the two together.

I’m writing this brief post because I found an R script from early 2012 in a random USB stick, in which coding superstar Jeff Hanson walked me through indexing, building loops, and functions in order to extract the species lineage ages for any phylogeny object.

The function below has a couple of loops that iterate through the mrca matrix and the branching times vector to produce a data frame of edge lengths for all the tree tips. The for loops make for slow processing if we have large trees but I was patient and that exact function created the main dataset that I analysed for a large part of my thesis. If you’re interested have a look at Jeff’s crisp programming to see how clever he was when coming up with the loops.

# extracting edge lengths from a phylo object

# define function

# note that the function requires ape to be installed

phylLages <- function(inputTree){

require(ape)

# numeric vector with the branching times for each tip

apedates<-branching.times(inputTree)

# matrix of tips and most recent common ancestors

mrca.matrix<-mrca(inputTree)

### preliminary procesing

Species = c()

Node_date = c()

### main processing

col_counter = 0

for (i in colnames(mrca.matrix)) {

  

  # keep track of which column we're in

  col_counter = col_counter + 1

  

  # print messege 

  cat(paste("Starting species ", as.character(col_counter), " out of ", as.character(length(colnames(mrca.matrix))), "\n", sep=""))

  

  # add species name to export vector

  Species = append(Species, i)

  

  # extract a vector of node numbers for the species we're up to

  currNodeVec = mrca.matrix[,col_counter]

  # working out what the node dates are for each node for the species we're working with

  currNodeDates = c()

  row_counter = 0

  for (currNodeNum in currNodeVec) {

    row_counter = row_counter + 1

    if (row_counter != col_counter) {

      currNodeDates = append(currNodeDates, apedates[as.character(currNodeNum)])

    }

  }

  # find out what the minimum number is

  minDate = min(currNodeDates)

  Node_date = append(Node_date, minDate)

  

        }

  exportFrame = data.frame(Species,Node_date)

  cat("DONE!\n")

  return(exportFrame)

}

A few years later I learned that the BioGeoBEARS R package includes the prt() function, which Nick Matzke wrote specifically to print the content of a tree into a tabular format, making all the implicit information explicit. This function is faster and provides much more information, and I was very relieved to see that the edge lengths it provides are the same Jeff and I obtained with the custom function.

Lets have a look. Both functions keep the tip order so comparing the results is easy. This examples uses the built in tree of bird families (Sibley and Ahlquist 1990) from the ape package.

plotted with ggtree

# comparing the edge lengths with BioGeoBEARS::prt

# load libraries

library(ape)

library(dplyr)

library(stringr)

library(BioGeoBEARS)

# phylogeny of bird families (comes with ape)

data("bird.families")

# get a data frame of edge lengths using the Jeff Hanson & Luis Verde function

jhEdges <- phylLages(bird.families)

# get a data frame of edge lengths using prt

birdPrintout <- prt(bird.families)

## extract the data for the tips only

# note: the tips in this particular tree are families not species 

birdPrintout_Tips <- birdPrintout %>% filter(str_detect(label,"Node")==FALSE) %>% 

                        select(Species=label,Node_date=edge.length)

# quick comparison

head(jhEdges) == head(birdPrintout_Tips)

 

Now we can do side-by-side glimpse of the two data frames, and time how long the two functions take to process the data.

> bind_cols(jhEdges,birdPrintout_Tips) %>% head

        Species Node_date       Species Node_date

1 Struthionidae      17.1 Struthionidae      17.1

2       Rheidae      17.1       Rheidae      17.1

3   Casuariidae       9.5   Casuariidae       9.5

4   Apterygidae       9.5   Apterygidae       9.5

5     Tinamidae      21.8     Tinamidae      21.8

6      Cracidae      19.8      Cracidae      19.8

system.time(phylLages(bird.families))

system.time(prt(bird.families))

prt is about 4.5x faster than the loops, and this becomes relevant as the number of tips increases or in the case of multiphylo objects with varying edge lengths so I recommend using it for any work on edge lengths.

==========================================

 

Regards,
Rui

> On Sep 11, 2017, at 11:50 PM, 张睿 <xiaor...@gmail.com> wrote:
>
> Hi Nick,
>
> Thank you for making a competent biogeographic tool BioGeoBEARS! It really helps for my phylogeographic analyses.
>
> Please add me to the google group and post the question I have here:
>
> How could I change the default colors for plotting per area?
>
> Ideally, I would like to something like:
>
> get_colors_for_numareas
> …...
>               area_colors = c(“my chosen color 1", "my chosen color 2", "my chosen color 3", "my chosen color 4", "my chosen color 5")
> ……
>
> Thank you very much!
>
> Regards,
> Rui
>