Modified version of Grafen's branch lengths

Last night - based on an R-sig-phylo request I wrote some code to compute a modified version of Grafen's edge lengths, in which the depth of each node is defined by the number of branching events above it - not the number of leaves.

However, this is not what I think the commenter wanted. Indeed, I believe he was looking for a transformation in which the height above the root of all nodes with a maximum path length to any tip of x was standardized for a given value of x.

Here is an example of this. First, our functions:

node.paths<-function(tree,node){
    d<-Descendants(tree,node,"children")
    paths<-as.list(d)
    while(any(d>Ntip(tree))){
        jj<-1
        new.paths<-list()
        for(i in 1:length(paths)){
            if(paths[[i]][length(paths[[i]])]<=Ntip(tree)){ 
                new.paths[[jj]]<-paths[[i]]
                jj<-jj+1
            } else {
                ch<-Descendants(tree,paths[[i]][length(paths[[i]])],
                    "children")
                for(j in 1:length(ch)){
                    new.paths[[jj]]<-c(paths[[i]],ch[j])
                    jj<-jj+1
                }
            }
        }
        paths<-new.paths
        d<-sapply(paths,function(x) x[length(x)])
    }
    paths
}
modified.Grafen<-function(tree,power=2){
    max.np<-function(tree,node){
        np<-node.paths(tree,node)
        if(length(np)>0) max(sapply(np,length)) else 0
    }
    nn<-1:(Ntip(tree)+tree$Nnode)
    h<-sapply(nn,max.np,tree=tree)+1
    h<-(h/max(h))^power
    edge.length<-vector()
    for(i in 1:nrow(tree$edge)) 
        edge.length[i]<-diff(h[tree$edge[i,2:1]])
    tree$edge.length<-edge.length
    tree
}

Now our trees:

library(phytools)
library(phangorn)
t1<-read.tree(text="((((A,B),C),E),(((F,G),H),I));")
t2<-read.tree(text="((((A,B),C,D),E),(((F,G),H),I));")
t3<-read.tree(text="((A,(B1,B2),C),(D,(E1,E2),F,(G1,G2)));")
## standard Grafen tree 1
plotTree(t1,lwd=1)

## modified Grafen
plotTree(modified.Grafen(t1),lwd=1)

## standard Grafen tree 2
plotTree(t2,lwd=1)

plotTree(modified.Grafen(t2),lwd=1)
nodelabels()

Note that nodes #11 & #14 have the same height because they have the same maximum path length to the tips.

## standard Grafen tree 3
plotTree(t3,lwd=1)

plotTree(modified.Grafen(t3,power=3),lwd=1)

Changing power just changes how the tree looks - not the fact that nodes a common maximum number of edges away from the present are at a common height.

Finally, here I show the plot with node labels giving the maximum depth (in nodes) of each internal node of the tree. We can see that with modified.Grafen nodes with the same index have the same height.

tree<-rtree(n=40)
tree$edge.length[sample(1:nrow(tree$edge),40)]<-0
tree<-di2multi(tree)
tree$edge.length<-NULL
max.np<-sapply(1:tree$Nnode+Ntip(tree),function(x,tree)
    max(sapply(node.paths(tree,x),length)),tree=tree)
plotTree(tree,lwd=1,fsize=0.7)
nodelabels(max.np,node=1:tree$Nnode+Ntip(tree),cex=0.7)

plotTree(modified.Grafen(tree,power=3),lwd=1,fsize=0.7)
nodelabels(max.np,node=1:tree$Nnode+Ntip(tree),cex=0.7)

That's all.

Cospeciation method in phytools

Today I have been working on a method to test the hypothesis that two trees are more topologically congruent than expected by chance. This is supposed to be a test for co-speciation.

This test has been used in the literature on a number of occasions (e.g., Logdon et al. 2011); however I'm unsure of the original reference as it doesn't seem to be an article by Rod Page, who has devised a number of methods to study cospeciation.

The method works by first calculating the distance between two input trees by a user specified criterion (e.g., Robinson-Foulds or SPR), then we simulate 100 or more stochastic, uncorrelated trees with labels matching our original phylogenies, finally the function computes a p-value for the test by comparison to this distribution. Note that I also have written plot and print methods for the function.

Here's an example. First, the functions:

cospeciation<-function(t1,t2,distance=c("RF","SPR"),
    method=c("simulation","permutation"),assoc=NULL,
    nsim=100,...){
    distance<-distance[1]
    if(!distance%in%c("RF","SPR")) distance<-"RF"
    method<-method[1]
    if(!method%in%c("simulation","permutation")) 
        method<-"simulation"
    if(is.null(assoc)){
        ## assume exact match
        tips<-intersect(t1$tip.label,t2$tip.label)
        assoc<-cbind(tips,tips)
    }
    if(any(!t1$tip.label%in%assoc[,1])) 
        t1<-drop.tip(t1,setdiff(t1$tip.label,assoc[,1]))
    if(any(!assoc[,1]%in%t1$tip.label))
        assoc<-assoc[assoc[,1]%in%t1$tip.label,]
    if(any(!t2$tip.label%in%assoc[,2])) 
        t2<-drop.tip(t2,setdiff(t2$tip.label,assoc[,2]))
    if(any(!assoc[,2]%in%t2$tip.label))
        assoc<-assoc[assoc[,2]%in%t2$tip.label,]    
    if(method=="permutation"){
        perm.labels<-function(tree){
            tree$tip.label<-sample(tree$tip.label)
            tree
        }
        tt1<-replicate(nsim,perm.labels(t1),simplify=FALSE)
        swap.t2<-t2
        swap.t2$tip.label<-sapply(t2$tip.label,function(x,y)
            y[which(y[,2]==x),1],y=assoc)
        tt2<-replicate(nsim,perm.labels(swap.t2),simplify=FALSE)
    } else {
        tt1<-pbtree(n=Ntip(t1),tip.label=t1$tip.label,
            nsim=nsim)
        swap.t2<-t2
        swap.t2$tip.label<-sapply(t2$tip.label,function(x,y)
            y[which(y[,2]==x),1],y=assoc)
        tt2<-pbtree(n=Ntip(t2),tip.label=swap.t2$tip.label,
            nsim=nsim)

    }
    if(distance=="SPR"){
        d.null<-mapply(SPR.dist,tt1,tt2)
        d<-SPR.dist(t1,swap.t2)
        P.val<-mean(c(d,d.null)<=d)
    } else {
        d.null<-mapply(RF.dist,tt1,tt2)
        d<-RF.dist(t1,swap.t2)
        P.val<-mean(c(d,d.null)<=d)
    }
    obj<-list(d=d,d.null=d.null,P.val=P.val,
        distance=if(distance=="SPR") "SPR" else "RF",
        method=if(method=="simulation") "simulation" else
        "permutation")
    class(obj)<-"cospeciation"
    obj
}

print.cospeciation<-function(x,...){
    cat(paste("\nCo-speciation test based on",x$distance,
        "distance.\n"))
    cat(paste("P-value obtained via",x$method,".\n\n"))
    if(x$distance=="SPR")
        cat(paste("   SPR distance:",x$d,"\n"))
    else
        cat(paste("   RF distance:",x$d,"\n"))
    cat(paste("   Mean(SD) from null: ",
        round(mean(x$d.null),1),"(",
        round(sd(x$d.null),1),")\n",sep=""))
    cat(paste("   P-value:",round(x$P.val,6),"\n\n"))
}

plot.cospeciation<-function(x,...){
    if(x$distance=="RF")
        p<-hist(x$d.null,breaks=seq(min(c(x$d,x$d.null))-3,
            max(c(x$d,x$d.null))+3,2),plot=FALSE)
    else if(x$distance=="SPR")
        p<-hist(x$d.null,breaks=seq(min(c(x$d,x$d.null))-1.5,
            max(c(x$d,x$d.null))+1.5,1),plot=FALSE)
    plot(p$mids,p$density,xlim=c(min(c(x$d,x$d.null))-2,
        max(c(x$d,x$d.null))+1),ylim=c(0,1.2*max(p$density)),
        type="n",
        xlab=paste(x$distance," distance (null by ",
        x$method,")",sep=""),ylab="relative frequency")
    y2<-rep(p$density,each=2)
    y2<-y2[-length(y2)]
    x2<-rep(p$breaks[2:length(p$breaks)-1],each=2)[-1]
    x3<-c(min(x2),x2,max(x2))
    y3<-c(0,y2,0)
    polygon(x3,y3,col=make.transparent("blue",0.3),
        border=FALSE)
    lines(p$breaks[2:length(p$breaks)-1],p$density,type="s")
    arrows(x$d,max(c(0.2*max(p$density),
        1.1*p$density[which(p$mids==x$d)])),
        x$d,0,lend="round",
        length=0.15,lwd=2,col="black")
    text(x$d-diff(par()$usr[1:2])/40,
        1.1*max(c(0.2*max(p$density),
        1.1*p$density[which(p$mids==x$d)])),
        "observed distance",
        srt=60,pos=4)
}

Now our trees. t1 & t2 were simulated with co-speciation and host-switching (via SPRs); whereas t3 is uncorrelated. Finally a1 maps the tips of t1 & t2; whereas t1 & t3 are related via a2.

library(phytools)
library(phangorn)
t1

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  R, F, Q, G, A, P, ...
## 
## Rooted; includes branch lengths.

t2

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  t.R, t.F, t.Q, t.G, t.A, t.P, ...
## 
## Rooted; includes branch lengths.

a1

##       tips tips 
##  [1,] "R"  "t.R"
##  [2,] "F"  "t.F"
##  [3,] "Q"  "t.Q"
##  [4,] "G"  "t.G"
##  [5,] "A"  "t.A"
##  [6,] "P"  "t.P"
##  [7,] "X"  "t.X"
##  [8,] "D"  "t.D"
##  [9,] "U"  "t.U"
## [10,] "M"  "t.M"
## [11,] "J"  "t.J"
## [12,] "H"  "t.H"
## [13,] "C"  "t.C"
## [14,] "N"  "t.N"
## [15,] "V"  "t.V"
## [16,] "E"  "t.E"
## [17,] "O"  "t.O"
## [18,] "S"  "t.S"
## [19,] "T"  "t.T"
## [20,] "Y"  "t.Y"
## [21,] "Z"  "t.Z"
## [22,] "I"  "t.I"
## [23,] "K"  "t.K"
## [24,] "W"  "t.W"
## [25,] "B"  "t.B"
## [26,] "L"  "t.L"

plot(cophylo(t1,t2,assoc=a1))

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

co1<-cospeciation(t1,t2,assoc=a1)
co1

## 
## Co-speciation test based on RF distance.
## P-value obtained via simulation .
## 
##    RF distance: 16 
##    Mean(SD) from null: 36.6(4.6)
##    P-value: 0.009901

plot(co1)

t3

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  r, w, o, g, n, t, ...
## 
## Rooted; includes branch lengths.

a2

##       LETTERS letters
##  [1,] "A"     "a"    
##  [2,] "B"     "b"    
##  [3,] "C"     "c"    
##  [4,] "D"     "d"    
##  [5,] "E"     "e"    
##  [6,] "F"     "f"    
##  [7,] "G"     "g"    
##  [8,] "H"     "h"    
##  [9,] "I"     "i"    
## [10,] "J"     "j"    
## [11,] "K"     "k"    
## [12,] "L"     "l"    
## [13,] "M"     "m"    
## [14,] "N"     "n"    
## [15,] "O"     "o"    
## [16,] "P"     "p"    
## [17,] "Q"     "q"    
## [18,] "R"     "r"    
## [19,] "S"     "s"    
## [20,] "T"     "t"    
## [21,] "U"     "u"    
## [22,] "V"     "v"    
## [23,] "W"     "w"    
## [24,] "X"     "x"    
## [25,] "Y"     "y"    
## [26,] "Z"     "z"

plot(cophylo(t1,t3,assoc=a2))

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

co2<-cospeciation(t1,t3,assoc=a2)
co2

## 
## Co-speciation test based on RF distance.
## P-value obtained via simulation .
## 
##    RF distance: 46 
##    Mean(SD) from null: 45.2(1.1)
##    P-value: 1

plot(co2)

OK, well, the plot method may still need some work, but you get the idea.