[R-sig-phylo] "Not of class phylo" error when lapply is used

Discussion:

jschenk

2018-04-09 14:46:37 UTC

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]
ages <- branching.times(phy)[as.character(NodeNumber)]
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")

#here is the same analysis conducted with multiple trees
multiTrees <- rmtree(20, 10)
class(multiTrees)
#I get the same error when I run my function across multiple trees
lapply(multiTrees, AgeDensity, species1="t3", species2="t15")

______________________________________________________________
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-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/

Liam J. Revell

2018-04-09 14:56:23 UTC

Permalink

Dear John.

You could try running .uncompressTipLabel on the "multiPhylo" object.
Let us know if that works.

All the best, Liam

Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston
& Profesor Asociado, Programa de Biología
Universidad del Rosario
web: http://faculty.umb.edu/liam.revell/

Post by jschenk
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]
ages <- branching.times(phy)[as.character(NodeNumber)]
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")
#here is the same analysis conducted with multiple trees
multiTrees <- rmtree(20, 10)
class(multiTrees)
#I get the same error when I run my function across multiple trees
lapply(multiTrees, AgeDensity, species1="t3", species2="t15")
______________________________________________________________
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]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/

jschenk

2018-04-10 02:37:53 UTC

Permalink

Thank you Graham, Liam, and Emmanuel for your suggestions. Applying the .uncompressTipLabel did the trick for my code. If anyone is interested in using the code to identify the 95% HPD for a single node from a set of chronograms, you can find the code here (https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R>) or below.

~John

library(ape)
library(coda)

#Function to estimate the age of a node given two members of a clade.
AgeDensity <- function(phy, species1, species2){
NodeNumber <- mrca(phy)[species1, species2]
ages <- branching.times(phy)[as.character(NodeNumber)]
return(as.numeric(ages))
}

#read in tree
tree <- read.nexus(file="")

#Apply function across the posterior distribution to obtain node ages. This will take some time to run, depending on the number of trees. It took me about 20 minutes.
NodeAges <- unlist(lapply(.uncompressTipLabel(tree), AgeDensity, "species1", "species2"))

#estimate the HPD interval for the node ages
HPD <- HPDinterval(as.mcmc(NodeAges), prob = 0.95)

#Density plot for 95% HPD
plot(density(as.numeric(NodeAges[which(NodeAges > HPD[1, 1] & NodeAges < HPD[1, 2])])), main="", las=1, xlab = "Million of years before present", lwd = 2)

______________________________________________________________
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

Post by Liam J. Revell
Dear John.
You could try running .uncompressTipLabel on the "multiPhylo" object. Let us know if that works.
All the best, Liam
Liam J. Revell, Associate Professor of Biology
University of Massachusetts Boston
& Profesor Asociado, Programa de Biología
Universidad del Rosario
web: http://faculty.umb.edu/liam.revell/

[[alternative HTML version deleted]]

_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/

Jacob Berv

2018-04-10 11:34:57 UTC

Permalink

This is cool — so this looks like it will work on a set of posterior trees even if they don’t all share the same topology?
Jake

Thank you Graham, Liam, and Emmanuel for your suggestions. Applying the .uncompressTipLabel did the trick for my code. If anyone is interested in using the code to identify the 95% HPD for a single node from a set of chronograms, you can find the code here (https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R><https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R>>) or below.
~John
library(ape)
library(coda)
#Function to estimate the age of a node given two members of a clade.
AgeDensity <- function(phy, species1, species2){
NodeNumber <- mrca(phy)[species1, species2]
ages <- branching.times(phy)[as.character(NodeNumber)]
return(as.numeric(ages))
}
#read in tree
tree <- read.nexus(file="")
#Apply function across the posterior distribution to obtain node ages. This will take some time to run, depending on the number of trees. It took me about 20 minutes.
NodeAges <- unlist(lapply(.uncompressTipLabel(tree), AgeDensity, "species1", "species2"))
#estimate the HPD interval for the node ages
HPD <- HPDinterval(as.mcmc(NodeAges), prob = 0.95)
#Density plot for 95% HPD
plot(density(as.numeric(NodeAges[which(NodeAges > HPD[1, 1] & NodeAges < HPD[1, 2])])), main="", las=1, xlab = "Million of years before present", lwd = 2)
______________________________________________________________
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 <http://sites.google.com/a/georgiasouthern.edu/schenk>
Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium <http://sites.google.com/a/georgiasouthern.edu/gasherbarium>

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>

jschenk

2018-04-10 12:49:45 UTC

Permalink

That’s right Jacob. Tracer will do the same thing if, and unfortunately only if, you predefine taxon sets in BEAUti. The problems with predefining clades in BEAUti are that (1) it is kind of a pain, (2) you often don’t know what the clades will be until after you conduct the analysis (as in my case where I am interested in shifts in diversification rates, which I won’t know until after I build the chronogram and conduct downstream analyses), and (3) Tracer doesn’t allow for any options on figure configuration.

~John

______________________________________________________________
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

Post by Jacob Berv
This is cool — so this looks like it will work on a set of posterior trees even if they don’t all share the same topology?
Jake

Thank you Graham, Liam, and Emmanuel for your suggestions. Applying the .uncompressTipLabel did the trick for my code. If anyone is interested in using the code to identify the 95% HPD for a single node from a set of chronograms, you can find the code here (https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R><https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R <https://github.com/johnjschenk/Rcode/blob/master/NodeAgeDensity.R>>) or below.
~John
library(ape)
library(coda)
#Function to estimate the age of a node given two members of a clade.
AgeDensity <- function(phy, species1, species2){
NodeNumber <- mrca(phy)[species1, species2]
ages <- branching.times(phy)[as.character(NodeNumber)]
return(as.numeric(ages))
}
#read in tree
tree <- read.nexus(file="")
#Apply function across the posterior distribution to obtain node ages. This will take some time to run, depending on the number of trees. It took me about 20 minutes.
NodeAges <- unlist(lapply(.uncompressTipLabel(tree), AgeDensity, "species1", "species2"))
#estimate the HPD interval for the node ages
HPD <- HPDinterval(as.mcmc(NodeAges), prob = 0.95)
#Density plot for 95% HPD
plot(density(as.numeric(NodeAges[which(NodeAges > HPD[1, 1] & NodeAges < HPD[1, 2])])), main="", las=1, xlab = "Million of years before present", lwd = 2)
______________________________________________________________
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 <http://sites.google.com/a/georgiasouthern.edu/schenk>
Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium <http://sites.google.com/a/georgiasouthern.edu/gasherbarium>

Post by jschenk
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]
ages <- branching.times(phy)[as.character(NodeNumber)]
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")
#here is the same analysis conducted with multiple trees
multiTrees <- rmtree(20, 10)
class(multiTrees)
#I get the same error when I run my function across multiple trees
lapply(multiTrees, AgeDensity, species1="t3", species2="t15")
______________________________________________________________
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 <http://sites.google.com/a/georgiasouthern.edu/schenk>
Herbarium website: sites.google.com/a/georgiasouthern.edu/gasherbarium <http://sites.google.com/a/georgiasouthern.edu/gasherbarium>
[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>

[[alternative HTML version deleted]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>

Emmanuel Paradis

2018-04-09 15:54:14 UTC

Permalink

Hi John,

See my comments below.

You may use:

NodeNumber <- getMRCA(phy, c(species1, species2))

this will be much more efficient for large trees.

Post by jschenk
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.

Post by jschenk
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)

Post by jschenk
#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.

Post by jschenk
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

Post by jschenk
______________________________________________________________
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]]
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Pour nous remonter une erreur de filtrage, veuillez vous rendre ici : http://f.security-mail.net/4074J0etF3M