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

2018-04-10 Thread jschenk
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 







> On Apr 10, 2018, at 7:34 AM, Jacob Berv  
> wrote:
> 
> 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
> 
>> On Apr 10, 2018, at 3:37 AM, jschenk > > wrote:
>> 
>> 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 
>> >  >) 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 
>>  
>> 
>> 
>> 
>> 
>> 
>> 
>> 
>>> On Apr 9, 2018, at 10:56 AM, Liam J. Revell >> > wrote:
>>> 
>>> 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/ 
>>> 
>>> 
>>> On 4/9/2018 9:46 AM, jschenk wrote:
 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 

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

2018-04-10 Thread 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

> On Apr 10, 2018, at 3:37 AM, jschenk  wrote:
> 
> 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 
>   >) 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 
>  
> 
> 
> 
> 
> 
> 
> 
>> On Apr 9, 2018, at 10:56 AM, Liam J. Revell  wrote:
>> 
>> 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/
>> 
>> On 4/9/2018 9:46 AM, jschenk wrote:
>>> 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 

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

2018-04-09 Thread jschenk
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 
) 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 







> On Apr 9, 2018, at 10:56 AM, Liam J. Revell  wrote:
> 
> 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/
> 
> On 4/9/2018 9:46 AM, jschenk wrote:
>> 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-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/


[[alternative HTML version deleted]]

___
R-sig-phylo 

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

2018-04-09 Thread Emmanuel Paradis

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/


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

2018-04-09 Thread 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/

On 4/9/2018 9:46 AM, jschenk wrote:

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-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/



___
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/