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"
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)
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")
Der Liam,
ReplyDeleteI 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!