Wednesday, May 4, 2016

as.phylo method for "simmap" and "multiSimmap" objects

I have just added some generic methods/functions as.phylo (method "simmap") and as.multiPhylo to strip the elements and class attributes from objects of class "simmap" and "multiSimmap". The reason for this, well, is because I discovered that the S3 method reorder for objects of class "simmap" was causing some problems for phangorn functions used internally in some of the new consensus tree functions of the phytools package.

This is really pretty simple. In the case of as.phylo, I just added the following function:

as.phylo.simmap<-function(x,...){
    x$maps<-NULL
    x$mapped.edge<-NULL
    if(!is.null(x$node.states)) x$node.states<-NULL
    if(!is.null(x$states)) x$states<-NULL
    if(!is.null(x$Q)) x$Q<-NULL
    if(!is.null(x$logL)) x$logL<-NULL
    if(!is.null(attr(x,"map.order"))) attr(x,"map.order")<-NULL
    class(x)<-setdiff(class(x),"simmap")
    x
}

and then exported an S3 method in NAMESPACE as follows:

S3method(as.phylo, simmap)

and that was all there was to it!

Here is a quick example using as.multiPhylo:

library(phytools)
packageVersion("phytools")
## [1] '0.5.29'
tree
## 
## Phylogenetic tree with 82 tips and 81 internal nodes.
## 
## Tip labels:
##  Anolis_ahli, Anolis_allogus, Anolis_rubribarbus, Anolis_imias, Anolis_sagrei, Anolis_bremeri, ...
## 
## Rooted; includes branch lengths.
head(ecomorph)
##        Anolis_ahli     Anolis_allogus Anolis_rubribarbus 
##                 TG                 TG                 TG 
##       Anolis_imias      Anolis_sagrei     Anolis_bremeri 
##                 TG                 TG                 TG 
## Levels: CG GB TC TG Tr Tw
map.trees<-make.simmap(tree,ecomorph,model="ER",nsim=10)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##             CG          GB          TC          TG          Tr          Tw
## CG -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145
## GB  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145  0.02314145
## TC  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145  0.02314145
## TG  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145  0.02314145
## Tr  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723  0.02314145
## Tw  0.02314145  0.02314145  0.02314145  0.02314145  0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##        CG        GB        TC        TG        Tr        Tw 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
map.trees
## 10 phylogenetic trees with mapped discrete characters
map.trees[[1]]
## 
## Phylogenetic tree with 82 tips and 81 internal nodes.
## 
## Tip labels:
##  Anolis_ahli, Anolis_allogus, Anolis_rubribarbus, Anolis_imias, Anolis_sagrei, Anolis_bremeri, ...
## 
## The tree includes a mapped, 6-state discrete character with states:
##  CG, GB, TC, TG, Tr, Tw
## 
## Rooted; includes branch lengths.
par(mfrow=c(4,3))
plot(map.trees,fsize=0.3,lwd=1)
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"
## no colors provided. using the following legend:
##        CG        GB        TC        TG        Tr        Tw 
##   "black"     "red"  "green3"    "blue"    "cyan" "magenta"

plot of chunk unnamed-chunk-3

trees<-as.multiPhylo(map.trees)
trees
## 10 phylogenetic trees
trees[[1]]
## 
## Phylogenetic tree with 82 tips and 81 internal nodes.
## 
## Tip labels:
##  Anolis_ahli, Anolis_allogus, Anolis_rubribarbus, Anolis_imias, Anolis_sagrei, Anolis_bremeri, ...
## 
## Rooted; includes branch lengths.
str(trees[[1]])
## List of 4
##  $ edge       : int [1:162, 1:2] 83 84 85 86 87 88 89 90 90 89 ...
##  $ Nnode      : int 81
##  $ tip.label  : chr [1:82] "Anolis_ahli" "Anolis_allogus" "Anolis_rubribarbus" "Anolis_imias" ...
##  $ edge.length: num [1:162] 0.268 0.209 0.402 0.826 0.768 ...
##  - attr(*, "class")= chr "phylo"
##  - attr(*, "order")= chr "cladewise"
par(mfrow=c(4,3))
nulo<-lapply(trees,plot,cex=0.3,no.margin=TRUE)

plot of chunk unnamed-chunk-4

Of course, the same could have been accomplished by simply setting the class attribute of the "simmap" objects to "phylo", "multiSimmap" to "multiPhylo", and so on. The only advantage of stripping away all the components associated with the object class is that these do not have to be continually copied & re-copied whenever the object is copied.

phytools can be updated thusly from GitHub using devtools:

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

Saturday, April 30, 2016

Plotting trees under a Gaussian / bell curve

I just submitted my Evolution 2016 talk title, and so I'd thought it might be fun to post the code that I used to generate the title slide image. The method is on the effect of consensus/average tree methods on PCMs, so I decided to (metaphorically) plot some trees under a Gaussian/bell curve. Just in case you should ever care to do the same, here is how I did it:

library(phytools)

mt<-phytools:::make.transparent

trees<-rmtree(n=5,N=100)

plot.new()

h<-sapply(trees,function(x) max(nodeHeights(x)))

x0<-10
x<-seq(0,20,by=0.01)
y<-100*exp(-(x-x0)^2/25)
par(mar=rep(0.1,4))
plot(x,y,type="l",axes=FALSE,xlab="",ylab="",col="transparent")

for(i in 1:length(trees)){
    x<-runif(n=1)*20
    y<-Inf
    while(y>100*exp(-((x+h[i]/2)-x0)^2/20)){
        x<-runif(n=1)*20
        y<-sample(1:90,1)
    }
    plotTree(trees[[i]],xlim=c(0,20)-x,ylim=c(0,100),
        tips=setNames(1:Ntip(trees[[i]])+y,
        trees[[i]]$tip.label),ftype="off",add=TRUE,lwd=1,
        color=mt("grey",0.4))
}

x0<-mean(par()$usr[1:2])
x<-seq(par()$usr[1],par()$usr[2],by=0.01)
y<-100*exp(-(x-x0)^2/25)
lines(x,y,lwd=2,col=mt("navy",0.3))

text("Consensus trees & their\nsuitability for macroevolutionary\nanalysis",
    x=mean(par()$usr[1:2]),y=1.3*mean(par()$usr[3:4]),cex=3.4)

text("Liam J. Revell (UMass-Boston)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4]),cex=1.8)
text("Klaus Schliep (UMass-Boston)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4])-3.6*strheight("X"),cex=1.8)
text("Cristián Hernández Ulloa(Univ. de Concepción)",x=mean(par()$usr[1:2]),
    y=0.6*mean(par()$usr[3:4])-7.2*strheight("X"),cex=1.8)

plot of chunk unnamed-chunk-1

If we want to export to PDF, we just have to flank the code above with the following lines, and to comment out plot.new(). That even allows us to use alternative fonts (although to do so, we first have to install the extrafonts library and take some other preliminary steps that can easily be found online).

pdf(file="Evolution-title.pdf",font="Calibri",width=10,height=7.5)
## ... our code
dev.off()

That's all there is to it.

Sunday, April 24, 2016

Update to cophylo to permit links in assoc for taxa not in the tree

I just added a small update (1, 2) to cophylo that permits links to be specified in assoc for taxa that are absent from the tree.

So, for instance:

library(phytools)
t1
## 
## Phylogenetic tree with 20 tips and 19 internal nodes.
## 
## Tip labels:
##  A, B, C, D, F, G, ...
## 
## Rooted; includes branch lengths.
t2
## 
## Phylogenetic tree with 23 tips and 22 internal nodes.
## 
## Tip labels:
##  t16, t3, t15, t2, t13, t6, ...
## 
## Rooted; includes branch lengths.
assoc
##       LETTERS      
##  [1,] "A"     "t17"
##  [2,] "B"     "t7" 
##  [3,] "C"     "t6" 
##  [4,] "D"     "t26"
##  [5,] "E"     "t21"
##  [6,] "F"     "t9" 
##  [7,] "G"     "t20"
##  [8,] "H"     "t23"
##  [9,] "I"     "t8" 
## [10,] "J"     "t18"
## [11,] "K"     "t12"
## [12,] "L"     "t10"
## [13,] "M"     "t11"
## [14,] "N"     "t1" 
## [15,] "O"     "t5" 
## [16,] "P"     "t4" 
## [17,] "Q"     "t19"
## [18,] "R"     "t24"
## [19,] "S"     "t22"
## [20,] "T"     "t13"
## [21,] "U"     "t2" 
## [22,] "V"     "t16"
## [23,] "W"     "t15"
## [24,] "X"     "t3" 
## [25,] "Y"     "t14"
## [26,] "Z"     "t25"

Both assoc[,1] has species not present in t1, and assoc[,2] has species not present in t2:

setdiff(assoc[,1],t1$tip.label)
## [1] "E" "I" "O" "R" "W" "Z"
setdiff(assoc[,2],t2$tip.label)
## [1] "t21" "t12" "t25"

Nonetheless (and this wouldn't have worked before):

obj<-cophylo(t1,t2,assoc)
## Rotating nodes to optimize matching...
## Done.
plot(obj)

plot of chunk unnamed-chunk-3

That's it.

The version of phytools with this update can be installed from GitHub as follows:

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