Method to compute consensus edge lengths (branch lengths for a consensus topology)

I just added a simple phytools function to compute consensus edge lengths - or, more accurately, edge lengths for a consensus topology - by various seemingly sensible methods, as follows:

Method 1: Compute the mean edge length for each edge in the consensus tree setting the length for each tree in which the edge is absent to zero. (Default setting. Function arguments method="mean.edge" and if.absent="zero".)

Method 2: Compute the mean edge length, but ignore trees in which the edge is absent. (Function arguments method="mean.edge" and if.absent="ignore".)

Method 3: Compute the non-negative least squares edge lengths on the consensus tree using the mean patristic distance matrix. (Function argument method="least.squares".) If the input trees are rooted & ultrametric, this can be used to produce a consensus tree that is also ultrametric.

By default the consensus tree method used is majority rule consensus (e.g., consensus(trees,p=0.5)); however the user can supply their own consensus tree using the optional argument, logically, consensus.tree.

Code for this new function can be viewed on phytools GitHub page, here.

The following is a quick demo using a real empirical data set in which the tree tip labels have been changed:

trees

## 512 phylogenetic trees

## method 1
t1<-consensus.edges(trees)
plotTree(t1,fsize=0.4)

## method 2
t2<-consensus.edges(trees,if.absent="ignore")
plotTree(t2,fsize=0.4)

## method 3
t3<-consensus.edges(trees,method="least.squares")
plotTree(t3,fsize=0.4)

## method 1, but with a 95% consensus tree
t1.p95<-consensus.edges(trees,consensus.tree=consensus(trees,p=0.95))
plotTree(t1.p95,fsize=0.4)

It's hard to compare the results from each method, but we could, for instance, correlate the patristic distance matrices that each consensus tree implies:

tips<-t1$tip.label
plot(cophenetic(t1)[tips,tips],cophenetic(t2)[tips,tips],
    xlab="Method 1",ylab="Method 2")
h<-2*max(nodeHeights(t1))
lines(c(0,h),c(0,h),lty="dashed",col="red",lwd=3)

plot(cophenetic(t1)[tips,tips],cophenetic(t3)[tips,tips],
    xlab="Method 1",ylab="Method 3")
lines(c(0,h),c(0,h),lty="dashed",col="red",lwd=3)

…. you get the idea.

In this case, at least, it seems to matter relatively little which approach is chosen; however this is an empirical result, so that should not be generally assumed.

Note that this was developed for rooted trees; but it may apply to unrooted trees. I haven't checked yet.

Not sure what other methods are out there for computing consensus edge lengths, nor how easy they would be to add. This is not exactly my area of expertise.

Simple function to match labels between trees

As part of something else that I'm working on I wanted to match labels between two trees in a manner analagous to the phytools function matchNodes (that is, a matrix with the corresponding node indices for the two trees, if nodes could be matched, and NA otherwise).

Here is what that very simple function looks like:

matchLabels<-function(tr1,tr2){
    foo<-function(x,y) if(length(obj<-which(y==x))>0) obj else NA
    M<-cbind(1:Ntip(tr1),sapply(tr1$tip.label,foo,y=tr2$tip.label))
    colnames(M)<-c("tr1","tr2")
    M
}

Here is how it works:

library(phytools)
t1<-rtree(n=14)
t2<-rtree(n=10)
par(mfrow=c(1,2))
plotTree(t1,offset=0.6)
tiplabels()
plotTree(t2,offset=0.6)
tiplabels()

matchLabels(t1,t2)

##     tr1 tr2
## t14   1  NA
## t8    2   3
## t1    3   4
## t6    4   7
## t11   5  NA
## t3    6   6
## t5    7   2
## t2    8  10
## t7    9   5
## t9   10   8
## t13  11  NA
## t10  12   9
## t4   13   1
## t12  14  NA

matchLabels(t2,t1)

##     tr1 tr2
## t4    1  13
## t5    2   7
## t8    3   2
## t1    4   3
## t7    5   9
## t3    6   6
## t6    7   4
## t9    8  10
## t10   9  12
## t2   10   8

Or, equivalently:

obj<-cophylo(t1,t2,rotate=FALSE)
plot(obj)
tiplabels.cophylo(which="left")
tiplabels.cophylo(which="right")

That's all there is to it.

New version of phytools (0.5-20) submitted to CRAN

I have just submitted a new version of phytools (0.5-20) to CRAN. Obviously, it may not be accepted immediately by the CRAN gatekeepers - but, if not, hopefully we can get it online within the next few days. In the meantime, this version can also be installed directly from GitHub as follows:

library(devtools) ## devtools must be installed
install_github("liamrevell/phytools")

Updates from the previous CRAN version of phytools include, but are perhaps not restricted to:

1. Update to read.newick to permit arguments to be passed by the user to scan.

2. An alternative density tree plotting function, now called densityTree (1, 2).

3. A new function (phylo.heatmap) for phylogenetic heat maps, suggested to me by a colleague Nate Swenson (1, 2, 3, 4, 5).

4. A new function (dotTree) for plotting dots of different sizes to represent a continuous character trait at the tips of the tree (1, 2, 3).

5. A new option to turn off the plotted points at the tips of the tree for objects of class "cophylo" (and other methods that use the same plotting function internally; 1, 2).

6. Update to the traitgram plotting function phenogram to permit it to be used with nodelabels, etc. from ape.

7. Update to phenogram to invisibly return the coordinates of the plotted tips.

8. A new function (rootedge.to.singleton) which converts a tree with a root edge to one with singleton nodes; and updates to nodeHeights and nodeheight to include the root edge length.

9. Some updates to the S3 plotting methods for objects of class "contMap" and "densityMap".

10. Some changes to the argument default values in reroot (here).

11. A version of dotTree for discretely valued characters (using different colors for different character states.

Here's a quick demo of the new phytools plotting function phylo.heatmap:

library(phytools)
packageVersion("phytools")

## [1] '0.5.20'

phylo.heatmap(tree,X,standardize=TRUE,pts=FALSE,lwd=2,
    labels=FALSE)

Comments welcome!

Bug fixes and updates for phytools; new name for density.tree to fix S3 method incompatibility

I'll be trying to get a new version of phytools on CRAN in advance of the workshop I will be teaching with Transmitting Science in Barcelona next week. Consequently, I just pushed a whole bunch of new updates to GitHub. Most of these pertain to minor bugs & fixes, but I also renamed the function density.tree as densityTree to rectify an unidentified (by me) S3 method incompatibility. This makes it even more confusing to distinguish this function from the similarly named, and similarly purposed, phangorn function densiTree. Both allow the user to visualize a posterior density of phylogenies from (say) a Bayesian analysis; however different plotting methods & options are used internally, and my function (now densityTree) can also plot mapped stochastic character histories. For example:

library(phytools)
packageVersion("phytools")

## [1] '0.5.19'

trees

## 100 phylogenetic trees

x

##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "b" "a" "b" "b" "b" "b" "b" "b" "a" "b" "b" "a" "a" "a" "a" "a" "b" "b" 
##   S   T   U   V   W   X   Y   Z 
## "a" "a" "a" "a" "a" "a" "a" "a"

## make stochastic character histories for x on trees
mtrees<-make.simmap(trees,x,message=FALSE)
mtrees

## 100 phylogenetic trees with mapped discrete characters

## plot posterior density using densityTree
colors<-setNames(c("red","blue"),sort(unique(x)))
densityTree(mtrees,colors=colors,nodes="inner",ftype="i",lwd=5,
    method="plotSimmap")
add.simmap.legend(colors=colors,prompt=FALSE,x=par()$usr[1]+0.1,
    y=par()$usr[4]-1)

After I submit to CRAN (hopefully soon) I will post a longer note documenting all the updates from the previous CRAN version. This version of phytools can be obtained from GitHub using the devtools package:

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