Co-phylogenetic plot for trees with no edge lengths

The phytools function plot.cophylo will no plot an object of class "cophylo" containing trees without edge lengths (upon user request).

Here's what I mean:

library(devtools)
install_github("liamrevell/phytools")

Now try it….

library(phytools)
t1<-rtree(n=26,tip.label=LETTERS,br=NULL)
t2<-rtree(n=26,tip.label=LETTERS,br=NULL)
obj<-cophylo(t1,t2)

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

plot(obj)

This was super easy. I just plugged the line:

if(is.null(tree$edge.length)) tree<-compute.brlen(tree)

into the internally used plotting function, phylogram, and that was all there was too it.

Small update to the small phytools function to fix malconformed trees

I just made a small update to the phytools function untangle, which basically tries to 'fix' malconformed "phylo" class objects that occasionally result from certain types of tree manipulation.

The update basically allows the returned "phylo" object to inherit the attributes of the input tree (except the ones that change in the fix) - which it otherwise may not depending on the method that is used to fix the malconformation.

Here's a quick demo:

library(phytools)
print(anole.trees,details=TRUE)

## 9 phylogenetic trees
## tree 1 : 18 tips
## tree 2 : 21 tips
## tree 3 : 11 tips
## tree 4 : 29 tips
## tree 5 : 26 tips
## tree 6 : 28 tips
## tree 7 : 17 tips
## tree 8 : 26 tips
## tree 9 : 19 tips

## estimate an MRP super tree
super.tree<-mrp.supertree(anole.trees)

## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"
## [1] "Best pscore so far: 174"

## The MRP supertree, optimized via pratchet(),
## has a parsimony score of 174 (minimum 174)

super.tree<-root(super.tree,outgroup="Iguana_iguana",resolve.root=TRUE) ## fails

## Warning in newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode: number of items
## to replace is not a multiple of replacement length

## Error in phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode: number of items to replace is not a multiple of replacement length

attributes(super.tree)

## $names
## [1] "edge"      "Nnode"     "tip.label"
## 
## $class
## [1] "phylo"
## 
## $order
## [1] "postorder"
## 
## $pscore
## [1] 174

super.tree<-untangle(super.tree,method="read.tree")
attributes(super.tree)

## $names
## [1] "edge"      "tip.label" "Nnode"    
## 
## $class
## [1] "phylo"
## 
## $order
## [1] "cladewise"
## 
## $pscore
## [1] 174

super.tree<-root(super.tree,outgroup="Iguana_iguana",resolve.root=TRUE) ## works
plotTree(super.tree,ftype="i",fsize=0.7)

That's it.

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!