Hi John,
See my comments below.
Le 09/04/2018 à 16:46, jschenk a écrit :
Hi Folks,
I have been banging my head against what appears to be an easy coding problem
for a while now and haven’t been able to hack my way out of it. I am running a
function to identify a posterior set of node ages for a particular node. The
function I wrote works just fine, but when I use lapply to sample across a
posterior distribution, I get the “not of class ‘phylo’” error, even when I
input a single tree of class phylo. I have tried every typical solution (e.g.,
reassigning class), but haven’t identified a solution. Does anyone know a
workaround?
Please see the example code below.
Thanks,
John
library(ape)
#simulate a single tree with 20 tips
simtree <- rtree(20)
#Make sure the tree exists
plot(simtree)
#Function I wrote to find a node and then tell me the age of the node. I
realize that the simulated tree is not ultrametric in this example, in real
life it will be - ultrametric trees also result in the same error
AgeDensity <- function(phy, species1, species2){
NodeNumber <- mrca(phy)[species1, species2]
You may use:
NodeNumber <- getMRCA(phy, c(species1, species2))
this will be much more efficient for large trees.
ages <- branching.times(phy)[as.character(NodeNumber)]
Beware that branching.times() is not meaningful for non-ultrametric
trees (as those generated by rtree). See this for instance:
R> tr <- rtree(2)
R> tr$edge.length
[1] 0.2098386 0.0353521
R> branching.times(tr)
3
0.0353521
R> branching.times(rotate(tr, 3))
3
0.2098386
If you really want to work with non-ultrametric trees, maybe you need to
use dist.nodes instead.
return(as.numeric(ages))
}
#check the class of the tree object, it will say that it is of class phylo
class(simtree)
#Run my function AgeDensity and it works just fine
AgeDensity(simtree, "t3", "t15")
#When I use the lapply function, I get an error that the object is not a of
class phylo, although I already verified that it is of class phylo.
lapply(simtree, AgeDensity, species1="t3", species2="t15")
You can build a list of trees simply:
simtree <- c(simtree)
where you can put several trees (as many as you want):
simtree <- c(rtree(20), rtree(20))
and add trees on an existing list:
simtree[[4]] <- rtree(20)
simtree[5:10] <- rmtree(6, 20)
#here is the same analysis conducted with multiple trees
multiTrees <- rmtree(20, 10)
Check-out the options of rmtree: what you did is similar to rmtree(N=20,
n=10) with N: number of trees, and n: number of tips.
class(multiTrees)
#I get the same error when I run my function across multiple trees
lapply(multiTrees, AgeDensity, species1="t3", species2="t15")
Indeed, but the error message is certainly different (and not only
because of the language set-up):
R> lapply(multiTrees, AgeDensity, species1="t3", species2="t15")
Error in mrca(phy)[species1, species2] : indice hors limites
Maybe you need to add a check in your AgeDensity function:
if (! species1 %in% phy$tip.label) stop("species1 not in the tree")
and the same for species2.
HTH
Best,
Emmanuel
______________________________________________________________
John J. Schenk, Ph.D.
Assistant Professor of Plant Biology
Georgia Southern University Herbarium (GAS), Curator
Department of Biology
4324 Old Register Road
Georgia Southern University
Statesboro, GA 30460-8042
Office: 2260 Biology Building
Office phone: (912) 478-0848
Lab website: sites.google.com/a/georgiasouthern.edu/schenk
Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Pour nous remonter une erreur de filtrage, veuillez vous rendre ici :
http://f.security-mail.net/4074J0etF3M
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
