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.