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")

1 comment:

  1. Der Liam,
    I am a big fan of phytools and have been exploring different approaches for ASR for my data dealing with niche specialization for ants. I have a question regarding a roadblock that I have been running into. I have selected a subset of BEAST trees from the posterior (n=1000) and my plan was to run a simmap using that random 1000 trees. I ran into several issues attempting to get that tree to run, but was able to run with the following:

    # Function to run make.simmap for a single tree
    run_make_simmap <- function(tree, trait_data) {
    make.simmap(tree, x = trait_data, model = "ER", Q = "mcmc")
    }

    # Initialize a list to store simulated mapped trees
    simmap_results <- vector("list", length = length(random.trees_pruned_for_mcmc))

    # Match binary trait data to tip labels for each tree in the multiphylo object
    for (i in seq_along(random.trees_pruned_for_mcmc)) {
    tree <- random.trees_pruned_for_mcmc[[i]]
    tip_labels <- tree$tip.label
    trait_data_named <- setNames(nesting_habit_cooperi_opport_mcmc, tip_labels)

    # Run make.simmap for the current tree
    simmap_results[[i]] <- run_make_simmap(tree, trait_data_named)
    }

    However, now I have a list and not a multiSimmap. Is there a way to merge that list into a single simmap and plot it using densityMap? Or I am missing a more straightforward approach for the original subset of trees? Thanks in advance for any feedback you could provide!

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.