Re-write of internally used function phyloDesign

I have totally re-written the internally used function, phyloDesign, which computes a design matrix for least-squares phylogeny estimation in the phytools function optim.phylo.ls, which does unweighted least-squares phylogeny estimation from a distance matrix.

Here's a simple demo using a DNA dataset from Jackman et al. (1999; Syst. Biol.).

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)

packageVersion("phytools")

## [1] '0.5.3'

library(phytools)
## load data
dna<-read.nexus.data(file="http://www.phytools.org/UMB2015/data/Jackman-etal.nex")
dna<-as.DNAbin(dna)
dna

## 55 DNA sequences in binary format stored in a list.
## 
## All sequences of same length: 1456 
## 
## Labels: Diplolaemus_darwinii Polychrus_acutirostris Anolis_acutus Anolis_cristatellus Anolis_krugi Anolis_stratulus ...
## 
## Base composition:
##     a     c     g     t 
## 0.339 0.261 0.118 0.281

## compute distance matrix
D<-dist.dna(dna)
## estimate tree, using NJ as a starting tree:
ls.fit<-optim.phylo.ls(D,phangorn::NJ(D))

## 1 set(s) of nearest neighbor interchanges. best Q so far = 0.187309346 
## 2 set(s) of nearest neighbor interchanges. best Q so far = 0.185881266 
## 3 set(s) of nearest neighbor interchanges. best Q so far = 0.1847435677 
## 4 set(s) of nearest neighbor interchanges. best Q so far = 0.1839579959 
## 5 set(s) of nearest neighbor interchanges. best Q so far = 0.1835163321 
## 6 set(s) of nearest neighbor interchanges. best Q so far = 0.1830769343 
## 7 set(s) of nearest neighbor interchanges. best Q so far = 0.1827496093 
## 8 set(s) of nearest neighbor interchanges. best Q so far = 0.18261176 
## 9 set(s) of nearest neighbor interchanges. best Q so far = 0.1823425143 
## 10 set(s) of nearest neighbor interchanges. best Q so far = 0.1821113405 
## 11 set(s) of nearest neighbor interchanges. best Q so far = 0.1819861917 
## 12 set(s) of nearest neighbor interchanges. best Q so far = 0.1819386781 
## 13 set(s) of nearest neighbor interchanges. best Q so far = 0.1819014877 
## 14 set(s) of nearest neighbor interchanges. best Q so far = 0.1814029139 
## 15 set(s) of nearest neighbor interchanges. best Q so far = 0.1813641191 
## 16 set(s) of nearest neighbor interchanges. best Q so far = 0.180996378 
## best Q score of 0.180996378 found after 16 nearest neighbor interchange(s).

plotTree(root(ls.fit,outgroup="Diplolaemus_darwinii",resolve.root=TRUE),
    fsize=0.8,ftype="i")

It's hard to summarize where the speed-up in phyloDesign (which is many-fold for reasonable sized trees) comes from. The original version was written with loops & logical statements. The new version with apply family functions and set comparisons. Basically, I'm a bit better at R programming nowadays than I was five years ago when I wrote the earlier version of this function. (But still not all that good….)

That's it!

Scale bar for cophylo plotting

I just added the option to the phytools function cophylo (for co-phylogenetic plotting) to show separate scale bars for the right & left plotted trees.

Here's a quick demo of how it works:

library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)

## Installing 1 packages: expm

## package 'expm' successfully unpacked and MD5 sums checked

library(phytools)

## Loading required package: ape
## Loading required package: maps
## 
##  # ATTENTION: maps v3.0 has an updated 'world' map.        #
##  # Many country borders and names have changed since 1990. #
##  # Type '?world' or 'news(package="maps")'. See README_v3. #

First, simulate some trees:

t1<-rtree(n=26,tip.label=LETTERS)
t2<-rtree(n=26,tip.label=LETTERS)

The basic plots, no scale bar, look as follows:

obj<-cophylo(t1,t2)

## Rotating nodes to optimize matching...
## Done.

plot(obj)

Now let's try with the same scale bar, 2.5 units in length, on either side:

plot(obj,scale.bar=rep(2.5,2))

We can use different scale bars:

plot(obj,scale.bar=c(2.5,3.5))

And, finally, if we set either to zero it turns off:

plot(obj,scale.bar=c(3,0))

Note that the latter plot might be misinterpreted as a single scale bar for both the right & left plotted trees - but this would be wrong (unless the two trees happen to have exactly the same total height) as the trees are always rescaled internally before plotting in this function!

That's it!

A few minor fixes & updates

I made a few minor updates to phytools in the past couple of days.

First, I fixed a bug in plotSimmap (and, consequently, in any function that uses plotSimmap internally, such as plotTree). The bug seems to result when the edges of the tree are neither in "cladewise" nor in "pruningwise" order, and I have fixed it by substituting the more standard "postorder" edge order (for post-order tree traversal) for ape's related "pruningwise" order. The fix can be seen here.

Second, I updated the phytools function collapse.to.star which collapses a user-specified subtree to a star (unsurprisingly) so that it now accepts trees without edge lengths. Details of this update can be seen here.

Finally, third, and perhaps most amusingly, I corrected the spelling of "neighber" to "neighbor" in a message printed to report the number of nearest-neighbor-interchanges required to find the best tree in the phytools function for unweighted least-squares phylogeny inference from a distance matrix, optim.phylo.ls.

That's it.

New, reasonably fast method to compute the patristic distance between a pair of species in a large tree

A R-sig-phylo question yesterday read as follows (very slightly paraphrased):

“Given a large newick tree e.g. from the new open tree of life (304959 tips) I want to calculate a pairwise distance between a small subset of species. The tricks I've seen so far [Ed. ape::cophenetic, which computes a N × N patristic distance for a tree with N tips] do not work on a tree of that size. I want to determine single species to single species distances without having an algorithm attempt to calculate the distance for all possible species (which is what ape is doing). Is there a package/trick to accomplish this?”

As a matter of fact, this can be done in a reasonably efficient way using the phytools function fastHeight, which computes the height above the root of the MRCA of any pair of species. That is because, for two species i and j, the patristic distane between them is simply the sum of the heights above the root for species i and j minus two times the height above the root of the common ancestor of i & j.

To show that this works fairly well, let's try to write a little custom function, fastDist, to compute this distance using the heights as mentioned:

library(phytools)
fastDist<-function(tree,sp1,sp2){
    fastHeight(tree,sp1,sp1)+fastHeight(tree,sp2,sp2)-
        2*fastHeight(tree,sp1,sp2)
}

Now, let's test it on a somewhat large tree:

## simulate tree at random
tree<-rtree(n=5000)
sp1<-sample(tree$tip.label)[1]
sp1

## [1] "t941"

sp2<-sample(tree$tip.label)[2]
sp2

## [1] "t2799"

## compute using cophenetic
system.time(d<-cophenetic(tree)[sp1,sp2])

##    user  system elapsed 
##    5.20    0.68    5.96

d

## [1] 11.34199

## compute using fastDist
system.time(d<-fastDist(tree,sp1,sp2))

##    user  system elapsed 
##    0.05    0.00    0.04

d

## [1] 11.34199

Having done that, let's go for broke with a tree of 304,959, just as described:

tree<-rtree(n=304959)
tree

## 
## Phylogenetic tree with 304959 tips and 304958 internal nodes.
## 
## Tip labels:
##  t51258, t83258, t60186, t263493, t117091, t210153, ...
## 
## Rooted; includes branch lengths.

sp1<-sample(tree$tip.label)[1]
sp1

## [1] "t51585"

sp2<-sample(tree$tip.label)[2]
sp2

## [1] "t244730"

## just to check that our MRCA is not the root
fastMRCA(tree,sp1,sp2)

## [1] 427604

## compute using fastDist
system.time(d<-fastDist(tree,sp1,sp2))

##    user  system elapsed 
##    4.80    1.34    6.21

d

## [1] 12.93951

Cool.

Just for the heck of it, I have also added this function to phytools & it is available from GitHub already.