Re: [R-sig-phylo] Estimating tip rates of molecular evolution

[Jacob Berv]

One way to do it with pGLS and root to tip path length may be to include
the number of nodes along a root to tip path as a control for node density
effects (maybe with some transformation). I don't know how well that will
work, but I am also experimenting with similar questions.

Another way is to take the tip rates from a relaxed clock estimate of
divergence times. CoEvol is much more powerful but as you note, not super
practical.

If you have a good number of species pairs, that may be the least
problematic way to do it-- it sounds like you have a large tree so maybe
that is optimal. there is an R package called sisters that may help with
it. https://rdrr.io/github/bomeara/sisters/

Best,
Jake Berv




On Thu, Nov 17, 2022, 5:38 AM Karla Shikev  wrote:

> Dear all,
>
> I want to test for a relationship between interspecific rates of evolution
> and a continuous predictor variable. I didn't want to use species pairs
> because they end up looking a lot of information. The method implemented in
> Colevol is very elegant, but is computationally unfeasible with large
> trees.
>
> One alternative could be to optimize BLs on a fixed topology and then use
> root-to-tip distances as a measure of rate of evolution, but I'm concerned
> about node density artifacts.
>
> Any suggestions would be greatly appreciated.
>
> Karla
>
> [[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 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] Data normalization

[Jacob Berv]

Interesting thought. I am not familiar with these analyses, but it would not 
surprise me if your intuition is correct.
J


> On Sep 19, 2021, at 1:00 PM, Ferenc Tibor Kagan  wrote:
> 
> Dear r-sig-phylo community
> 
> I am writing to you in hopes of you giving me your inputs on a specific topic.
> 
> I noticed a rise of use of PCMs when it comes to gene expression data lately. 
> Many of these studies before fitting a specific model to their expression 
> data do several normalization steps. The common steps in order are to 
> normalize for sequencing depth and gene length, normalize in between 
> replicates within species and finally to normalize across species. For within 
> species and across species I have seen TMM normalization method being used 
> (from edgeR package) or batch effect removal (f.ex. Combat-seq function from 
> sva package).
> 
> My concern is the final normalization step, namely to normalize continuous 
> data across species before model fitting. By doing so wouldn't one minimize 
> the phylogenetic signal present in the dataset, therefore affecting the best 
> fitting model?
> 
> Thank you in advance for your answer.
> 
> Best regards,
> Ferenc Kagan
> 
>   [[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/

Re: [R-sig-phylo] Ancestral state reconstruction, polytomies and the absence of branch lengths

[Jacob Berv]

Hi Diego,

If your tree has polytomies and the branches have no length information, I am 
not sure if likelihood model-based reconstruction is right.

The ace function (and all other similar functions based on likelihood) will 
assume that your branch lengths represent some kind biological information 
(usually time), and the degree to which characters evolve will be related to 
the branch length. I therefore might prefer parsimony in this situation to get 
a general sense of what is going on.

As an experiment, you could assign all branches to be equal (e.g. = 1), but 
then you are imposing this assumption on the analysis. In order to evaluate the 
effects of this assumption, and if you are specifically interested in 
reconstructing the ancestral states for a particular node, you could generate a 
few sets of simulated branch lengths (which you could then pass to ace), to see 
if the reconstruction is robust to a range of variation in branch lengths. You 
could also replicate using multi2di to see how robust your reconstruction is to 
alternative resolutions.

Best,
Jake


> On Sep 15, 2021, at 2:16 AM, Diego Almeida-Silva 
>  wrote:
> 
> Dear colleagues,
> 
> I am trying to estimate ancestral states for two binary characters in a
> full morphology-based phylogenetic tree. I am using ace function,
> type="discrete" and kappa=0 as arguments. Unfortunately, I am dealing with
> two problems about this issue: my topology have three polytomies and tree
> branches have no lengths.
> 
> To perform ancestral state estimation, I firstly used the multi2di
> function. However, this procedure creates artificial nodes in the topology.
> This becomes a new problem, since states at these new nodes are also being
> estimated.
> 
> Is there another way to deal with this situation? There are functions (or
> packages) other than those of phytools that are able to deal with
> polytomies and the absence of branch lengths in ancestral state
> reconstruction?
> 
> Thanks in advance!
> Diego
> 
>   [[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/

[R-sig-phylo] identifying phylogenetically equivalent nodes

[Jacob Berv]

Dear R-sig-phylo,

Over the weekend, I asked Liam Revell if he had a solution to use matchNodes 
for a particular problem I’m trying to solve—finding all phylogenetically 
equivalent nodes when comparing trees that have uneven taxon samples and 
different topologies. Liam was kind enough to take some time to write a blog 
post about this, and got me started with some code

http://blog.phytools.org/2021/02/on-matching-nodes-between-trees-using.html

On it’s face this seems like a simple problem, but I’m running into some issues 
and thought I would reach out to the broader group. The code linked above seems 
to work, but only for comparing trees that start out as topologically 
identical. For my purposes, I’m trying to match nodes from a given a reference, 
to nodes in and across several hundred gene trees that differ in topology and 
taxon sample relative to the reference.

Here is a function definition based on Liam’s example

#function to match nodes from consensus 
#to individual gene trees with uneven sampling
#derived from Liam Revell's example-- need to 
testmatch_phylo_nodes<-function(t1, t2){
  ## step one drop tips
  t1p<-drop.tip(t1,setdiff(t1$tip.label, t2$tip.label))
  t2p<-drop.tip(t2,setdiff(t2 $tip.label, t1$tip.label))
  
  ## step two match nodes "descendants"
  M<-matchNodes(t1p,t2p)
  
  ## step two match nodes "distances"
  M1<-matchNodes(t1,t1p,"distances")
  M2<-matchNodes(t2,t2p,"distances")
  
  ## final step, reconcile
  MM<-matrix(NA,t1$Nnode,2,dimnames=list(NULL,c("left","right")))
  
  for(i in 1:nrow(MM)){
MM[i,1]<-M1[i,1]
nn<-M[which(M[,1]==M1[i,2]),2]
if(length(nn)>0){   
MM[i,2]<-M2[which(M2[,2]==nn),1]
}   
  }
  return(MM)
}


When t1 and t2 are trees that have topological conflicts, this function returns 
an error: 

Error in MM[i, 2] <- M2[which(M2[, 2] == nn), 1] : 
  replacement has length zero

I think(?) this happens because a particular node doesn’t exist in one or the 
other trees, and it returns integer(0) at that line — but I’m not sure I really 
understand what is going on here.


I modified Liam’s code slightly to get it to run without error in the described 
case, by making it conditional on that particular line:


Modified version

#function to match nodes from consensus 
#to individual gene trees with uneven sampling
#derived from Liam Revell's example-- need to test
match_phylo_nodes<-function(t1, t2){
## step one drop tips
t1p<-drop.tip(t1,setdiff(t1$tip.label, t2$tip.label))
t2p<-drop.tip(t2,setdiff(t2 $tip.label, t1$tip.label))

## step two match nodes "descendants"
M<-matchNodes(t1p,t2p)

## step two match nodes "distances"
M1<-matchNodes(t1,t1p,"distances")
M2<-matchNodes(t2,t2p,"distances")

## final step, reconcile
MM<-matrix(NA,t1$Nnode,2,dimnames=list(NULL,c("left","right")))

for(i in 1:nrow(MM)){
MM[i,1]<-M1[i,1]
nn<-M[which(M[,1]==M1[i,2]),2]
if(length(nn)>0){   
if(length(which(M2[,2]==nn))>0){
MM[i,2]<-M2[which(M2[,2]==nn),1]
}
} else {
}   
}
return(MM)  
}


I’ve been experimenting with this and some downstream code for the last few 
days, but I’ve run into some weird inconsistent results (not easily summarized) 
that make me think that this function is not working as intended.

I was wondering — have any of you dealt with a similar problem? In principle 
this seems like it should be similar to concordance analysis, but I care less 
about identifying the proportion of nodes that exist in gene trees given a 
reference, and instead I need the actual node numbers in a given gene tree that 
are phylogenetically equivalent to particular nodes in a reference. Happy to 
try to hack away at something… 


Best,
Jake Berv





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

Re: [R-sig-phylo] model averaging for discrete character evolution

[Jacob Berv]

Accidentally replied directly to Marguerite, instead of to the list: 

> On Aug 9, 2019, at 10:01 AM, Jacob Berv  
> wrote:
> 
> Great points Marguerite. I agree that the idea of using ‘all possible’ models 
> is probably a bit too overzealous— and as Brian pointed out, some models will 
> probably generate totally bogus results. For example, if you had three states 
> and only allowed one transition rate between two states-- clearly that would 
> not make any sense.
> 
> But perhaps it might be useful to explore a reasonable set of models which 
> have plausible biological interpretations — like the ARD, ER, SYM models 
> typically available (and perhaps a few custom models).
> 
> I was thinking that model averaging in this context may be a useful way to 
> integrate out some of the uncertainty in ‘model choice.’ AIC/BIC ranking is 
> one way to do that, but in my mind ‘best statistical fit' may not be 
> synonymous with ‘best,' at least without some hypothesis of the underlying 
> biology. Though I suppose if you are weighting by AIC/BIC then you are still 
> looking at 'statistical fit,’ just in another way.
> 
> J
> 
> 
>> On Aug 9, 2019, at 2:06 AM, Marguerite Butler > <mailto:mbutler...@gmail.com>> wrote:
>> 
>> Aloha Jacob, Liam, Brian,
>> 
>> To add to Brian's and Liam's excellent comments - just remember that any 
>> type of ancestral state reconstruction of any form is just a weighted 
>> average of the information that you put into the problem. If you have no 
>> fossil data or data back in time, it's a weighted average of the tips 
>> dictated by the model and/or the algorithm that you use.  So point one is 
>> keep that in mind. 
>> 
>> Second, along the lines of Brian's comments, it is always worth looking at 
>> the model outputs - first of all, are the estimates different?  If not, then 
>> there is no reason to model average, the solution is robust to the model.  
>> If they are different, then would it make more sense to interpret each model 
>> results one by one and think about what the scenarios imply? It might be 
>> useful to consider the most distinct alternative hypotheses, for example.   
>> Of course I always only include in the model set the ones that have a 
>> biological interpretation (the most robust tests have the most distinct but 
>> biologically interpretable hypotheses). It is a way of bracketing the 
>> interpretation of an answer.  
>> 
>> Often model averaging will result in parameter estimates that are 
>> intermediate in value - but if none of the models returned that estimate, 
>> then what is it exactly? It's a point that is best fit for nothing.  But the 
>> solution along an entire tree is a solution that is part of a set (a single 
>> evolutionary scenario).  Mixing and matching from different solutions may 
>> not make sense. 
>> 
>> So I've personally never found a situation where model averaging is useful 
>> and sometimes can lead to non-critical thinking.  But if you find something, 
>> I'm all ears!
>> 
>> Best of luck,
>> Marguerite
>> 
>> On Thu, Aug 8, 2019 at 10:14 PM Jacob Berv > <mailto:jakeberv.r.sig.ph...@gmail.com>> wrote:
>> Cool- thanks for the great comments Liam and Brian, as always!
>> J
>> 
>> > On Aug 7, 2019, at 11:55 AM, Liam Revell > > <mailto:liam.rev...@umb.edu>> wrote:
>> > 
>> > Hi Jake.
>> > 
>> > [Edit: I see just now that Brian has also responded to this inquiry. I 
>> > have no doubt that his message is more insightful than mine - but I'll 
>> > nonetheless send what I was writing anyway just in case it contributes 
>> > anything useful to the discussion.]
>> > 
>> > If you're simply interested in the states at nodes, you might consider 
>> > just multiplying the marginal reconstructions under maximum likelihood 
>> > by the Akaike weights of each model & summing them.
>> > 
>> > The former are probabilities that the node is in each state conditioned 
>> > on a model, and the latter are the probabilities that each of the models 
>> > is the best of the set. If the models in the set genuinely comprise all 
>> > possible ways in which your character could have evolved (they don't, 
>> > but still), then the model weighted average marginal reconstructions 
>> > should give the total probability that each node is in each state.
>> > 
>> > If you really want to do stochastic mapping, you might consider some 
>> > kind of rjMCMC in which you sample models & transition rates

Re: [R-sig-phylo] model averaging for discrete character evolution

[Jacob Berv]

Cool- thanks for the great comments Liam and Brian, as always!
J

> On Aug 7, 2019, at 11:55 AM, Liam Revell  wrote:
> 
> Hi Jake.
> 
> [Edit: I see just now that Brian has also responded to this inquiry. I 
> have no doubt that his message is more insightful than mine - but I'll 
> nonetheless send what I was writing anyway just in case it contributes 
> anything useful to the discussion.]
> 
> If you're simply interested in the states at nodes, you might consider 
> just multiplying the marginal reconstructions under maximum likelihood 
> by the Akaike weights of each model & summing them.
> 
> The former are probabilities that the node is in each state conditioned 
> on a model, and the latter are the probabilities that each of the models 
> is the best of the set. If the models in the set genuinely comprise all 
> possible ways in which your character could have evolved (they don't, 
> but still), then the model weighted average marginal reconstructions 
> should give the total probability that each node is in each state.
> 
> If you really want to do stochastic mapping, you might consider some 
> kind of rjMCMC in which you sample models & transition rates from their 
> joint posterior distribution and then generate stochastic character maps 
> based on this sample. This is not implemented in phytools::make.simmap 
> (it does MCMC, but only given a model), but is not to hard to envision 
> doing, so long as you are careful about designing the rjMCMC.
> 
> Now to read Brian's answer
> 
> All the best, Liam
> 
> Liam J. Revell
> Associate Professor, University of Massachusetts Boston
> Profesor Asistente, Universidad Católica de la Ssma Concepción
> web: http://faculty.umb.edu/liam.revell/ 
> <http://faculty.umb.edu/liam.revell/>, http://www.phytools.org 
> <http://www.phytools.org/>
> 
> Academic Director UMass Boston Chile Abroad (starting 2019):
> https://www.umb.edu/academics/caps/international/biology_chile 
> <https://www.umb.edu/academics/caps/international/biology_chile>
> 
> On 8/7/2019 11:32 AM, Jacob Berv wrote:
>> Dear R-sig phylo
>> 
>> I’ve been running a few discrete character Mk type models using phytool's 
>> SIMMAP — and I had the idea that it might be useful to try model averaging 
>> across posterior probabilities for node states.
>> 
>> Might this make sense to do, to avoid problems associated with model ranking 
>> via AIC? Ie, average the node state probabilities based on AIC weights? Is 
>> there some fundamental problem with this?
>> 
>> I could imagine generating some code to generate all possible transition 
>> models given a set of N states, and then rather than ranking with AIC, model 
>> averaging for parameter estimates (though now that I think about it, not 
>> sure how one might reasonably average a symmetrical rate and an asymmetrical 
>> rate).
>> 
>> Best,
>> Jake
>> 
>>  [[alternative HTML version deleted]]
>> 
>> ___
>> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
>> <mailto:R-sig-phylo@r-project.org>
>> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylodata=02%7C01%7Cliam.revell%40umb.edu%7Ce85015cb8c314c684bdc08d71b4c8fee%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C637007887976647805sdata=F8qdVoT17J%2BuC6wgvFktdI1id%2Bf37hDLe1W17z6OhTw%3Dreserved=0
>>  
>> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylodata=02%7C01%7Cliam.revell%40umb.edu%7Ce85015cb8c314c684bdc08d71b4c8fee%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C637007887976647805sdata=F8qdVoT17J%2BuC6wgvFktdI1id%2Bf37hDLe1W17z6OhTw%3Dreserved=0>
>> Searchable archive at 
>> https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2Fdata=02%7C01%7Cliam.revell%40umb.edu%7Ce85015cb8c314c684bdc08d71b4c8fee%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C637007887976657803sdata=ogCaep2DN92bk6%2FZwOOcpmHaCPQwX7BzNgKo2UMcWtw%3Dreserved=0
>>  
>> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2Fdata=02%7C01%7Cliam.revell%40umb.edu%7Ce85015cb8c314c684bdc08d71b4c8fee%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C637007887976657803sdata=ogCaep2DN92bk6%2FZwOOcpmHaCPQwX7BzNgKo2UMcWtw%3Dreserved=0>
>> 


[[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] model averaging for discrete character evolution

[Jacob Berv]

Dear R-sig phylo

I’ve been running a few discrete character Mk type models using phytool's 
SIMMAP — and I had the idea that it might be useful to try model averaging 
across posterior probabilities for node states.

Might this make sense to do, to avoid problems associated with model ranking 
via AIC? Ie, average the node state probabilities based on AIC weights? Is 
there some fundamental problem with this?

I could imagine generating some code to generate all possible transition models 
given a set of N states, and then rather than ranking with AIC, model averaging 
for parameter estimates (though now that I think about it, not sure how one 
might reasonably average a symmetrical rate and an asymmetrical rate).

Best,
Jake

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

Re: [R-sig-phylo] aov.phylo vs phylANOVA

[Jacob Berv]

Sorry— I meant aov.phylo (geiger function), not phly.aov. 
J

> On Nov 15, 2018, at 3:00 PM, Jacob Berv  
> wrote:
> 
> Hi everyone — 
> It seems like the solution is that I was looking at the wrong p-value, as is 
> often the case!
> 
> From Liam:
> "The issue may be that you need to look at the p-value from the internal 
> print-out of aov.phylo not from the summary."
> 
> I was using the summary() function to examine the model output from geiger. 
> When applying this function to the phyl.aov output object, the p-values 
> appear to be from a regular anova, not the phylogenetic anova. To access the 
> phylogenetic anova p value, you have to look inside the object generated by 
> phyl.aov:
> 
> attr(obj, ’summary’)
> 
> This will display the appropriate "Pr(>F) given phy"
> 
> Jake
> 
> 
>> On Nov 15, 2018, at 12:53 PM, Liam J. Revell  wrote:
>> 
>> Hi Jacob.
>> 
>> As far as I know, aov.phylo and phylANOVA should be doing more or less the 
>> same thing. With random data if I run enough simulations for the null 
>> distribution of F the P-values of the two different implementations come out 
>> almost exactly the same. One difference that I noted is that if you give 
>> either method data vectors without taxon names both work, but only phylANOVA 
>> gives a warning.
>> 
>> Please send along the data that has generated this incongruency if you are 
>> unable to figure it out.
>> 
>> All the best, Liam
>> 
>> Liam J. Revell
>> Associate Professor, University of Massachusetts Boston
>> Profesor Asistente, Universidad Católica de la Ssma Concepción
>> web: http://faculty.umb.edu/liam.revell/, http://www.phytools.org
>> 
>> On 11/15/2018 2:30 PM, Jacob Berv wrote:
>>> Dear R-sig-phylo,
>>> I was wondering if anyone on here might be able to help me understand the 
>>> difference between phytool’s implementation of phylogenetic ANOVA and 
>>> geiger’s implementation. From the respective documentation, it seems that 
>>> both approaches rely on and cite the same reference:
>>> Garland T Jr, AW Dickerman, CM Janis, and JA Jones. 1993. Phylogenetic 
>>> analysis of covariance by computer simulation. Systematic Biology 
>>> 42(3):265-292.
>>> Both seem to have a similar approach, at least as it is described in their 
>>> respective documentations, and both seem to rely on character simulations 
>>> to derive their p values. It seems aov.phylo uses sim.char() and phylANOVA 
>>> uses fastBM() for their simulations internally.
>>> On Liam’s blog, he indicates that these tests are the same, except that 
>>> phylANOVA additionally performs post-hoc tests.
>>> http://blog.phytools.org/2013/02/updated-phylanova.html
>>> However, running some of my data through both of these tests is generating 
>>> totally different results (aov.phylo detecting significant differences 
>>> where phylANOVA does not, with p values differing by 5 orders of magnitude.
>>> Running my same test data~group through a pgls also generates a result 
>>> comparable to what I get from phylANOVA — so it seems like perhaps 
>>> aov.phylo is the outlier?
>>> Best,
>>> Jake Berv
>>> ___
>>> 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/

Re: [R-sig-phylo] aov.phylo vs phylANOVA

[Jacob Berv]

Hi everyone — 
It seems like the solution is that I was looking at the wrong p-value, as is 
often the case!

>From Liam:
"The issue may be that you need to look at the p-value from the internal 
print-out of aov.phylo not from the summary."

I was using the summary() function to examine the model output from geiger. 
When applying this function to the phyl.aov output object, the p-values appear 
to be from a regular anova, not the phylogenetic anova. To access the 
phylogenetic anova p value, you have to look inside the object generated by 
phyl.aov:

attr(obj, ’summary’)

This will display the appropriate "Pr(>F) given phy"

Jake


> On Nov 15, 2018, at 12:53 PM, Liam J. Revell  wrote:
> 
> Hi Jacob.
> 
> As far as I know, aov.phylo and phylANOVA should be doing more or less the 
> same thing. With random data if I run enough simulations for the null 
> distribution of F the P-values of the two different implementations come out 
> almost exactly the same. One difference that I noted is that if you give 
> either method data vectors without taxon names both work, but only phylANOVA 
> gives a warning.
> 
> Please send along the data that has generated this incongruency if you are 
> unable to figure it out.
> 
> All the best, Liam
> 
> Liam J. Revell
> Associate Professor, University of Massachusetts Boston
> Profesor Asistente, Universidad Católica de la Ssma Concepción
> web: http://faculty.umb.edu/liam.revell/, http://www.phytools.org
> 
> On 11/15/2018 2:30 PM, Jacob Berv wrote:
>> Dear R-sig-phylo,
>> I was wondering if anyone on here might be able to help me understand the 
>> difference between phytool’s implementation of phylogenetic ANOVA and 
>> geiger’s implementation. From the respective documentation, it seems that 
>> both approaches rely on and cite the same reference:
>> Garland T Jr, AW Dickerman, CM Janis, and JA Jones. 1993. Phylogenetic 
>> analysis of covariance by computer simulation. Systematic Biology 
>> 42(3):265-292.
>> Both seem to have a similar approach, at least as it is described in their 
>> respective documentations, and both seem to rely on character simulations to 
>> derive their p values. It seems aov.phylo uses sim.char() and phylANOVA uses 
>> fastBM() for their simulations internally.
>> On Liam’s blog, he indicates that these tests are the same, except that 
>> phylANOVA additionally performs post-hoc tests.
>> http://blog.phytools.org/2013/02/updated-phylanova.html
>> However, running some of my data through both of these tests is generating 
>> totally different results (aov.phylo detecting significant differences where 
>> phylANOVA does not, with p values differing by 5 orders of magnitude.
>> Running my same test data~group through a pgls also generates a result 
>> comparable to what I get from phylANOVA — so it seems like perhaps aov.phylo 
>> is the outlier?
>> Best,
>> Jake Berv
>> ___
>> 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/

[R-sig-phylo] aov.phylo vs phylANOVA

[Jacob Berv]

Dear R-sig-phylo,
I was wondering if anyone on here might be able to help me understand the 
difference between phytool’s implementation of phylogenetic ANOVA and geiger’s 
implementation. From the respective documentation, it seems that both 
approaches rely on and cite the same reference:

Garland T Jr, AW Dickerman, CM Janis, and JA Jones. 1993. Phylogenetic analysis 
of covariance by computer simulation. Systematic Biology 42(3):265-292.

Both seem to have a similar approach, at least as it is described in their 
respective documentations, and both seem to rely on character simulations to 
derive their p values. It seems aov.phylo uses sim.char() and phylANOVA uses 
fastBM() for their simulations internally. 

On Liam’s blog, he indicates that these tests are the same, except that 
phylANOVA additionally performs post-hoc tests. 
http://blog.phytools.org/2013/02/updated-phylanova.html

However, running some of my data through both of these tests is generating 
totally different results (aov.phylo detecting significant differences where 
phylANOVA does not, with p values differing by 5 orders of magnitude. 

Running my same test data~group through a pgls also generates a result 
comparable to what I get from phylANOVA — so it seems like perhaps aov.phylo is 
the outlier?

Best,
Jake Berv
___
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] plot heatmap at nodes

[Jacob Berv]

Ah, not really. I am thinking more of an actual matrix, where matrix entries 
represent support values, and are colored proportionally. Ie - something like:

image(t(matrix(data=c(95, 60,75,85), nrow=2, ncol=2)), col=heat.colors(4), 
yaxt='n', xaxt='n')

Jake


> On Aug 29, 2018, at 8:06 PM, Liam J. Revell  wrote:
> 
> Hi Jake. Can you draw a clearer picture of what you'd like to do? Would it 
> look something like the error bars that phytools can compute for contMap 
> object? (e.g., 
> http://blog.phytools.org/2017/02/function-to-add-error-bars-to-contmap.html)
> 
> All the best, Liam
> 
> Liam J. Revell
> Associate Professor, University of Massachusetts Boston
> Profesor Asistente, Universidad Católica de la Ssma Concepción
> web: http://faculty.umb.edu/liam.revell/, http://www.phytools.org
> 
> On 8/29/2018 7:12 PM, Jacob Berv wrote:
>> Dear R-sig-phylo,
>> Does there exist a function to plot heatmaps at internal nodes of a tree? 
>> For example, to illustrate support values from different phylogenetic 
>> analyses — ie given a set of 3x3 matricies of support values for particular 
>> nodes.
>> It seems like there may be a partial (and definitely workable) solution 
>> here, using ggtree— 
>> https://rgriff23.github.io/2017/05/11/primate-phylogeny-ggtree.html 
>> <https://rgriff23.github.io/2017/05/11/primate-phylogeny-ggtree.html>, if 
>> such heatmaps are saved as raster images (about 1/2 way down). This may be 
>> the easiest solution, but it would be preferable to keep everything in 
>> vector format and all in R.
>> I guess I could do this by calling image() for the individual heatmaps and 
>> then just linking them up to nodes of interest in illustrator (but that 
>> feels like cheating).
>> Jake
>>  [[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/

[R-sig-phylo] plot heatmap at nodes

[Jacob Berv]

Dear R-sig-phylo,

Does there exist a function to plot heatmaps at internal nodes of a tree? For 
example, to illustrate support values from different phylogenetic analyses — ie 
given a set of 3x3 matricies of support values for particular nodes.

It seems like there may be a partial (and definitely workable) solution here, 
using ggtree— 
https://rgriff23.github.io/2017/05/11/primate-phylogeny-ggtree.html 
, if such 
heatmaps are saved as raster images (about 1/2 way down). This may be the 
easiest solution, but it would be preferable to keep everything in vector 
format and all in R.

I guess I could do this by calling image() for the individual heatmaps and then 
just linking them up to nodes of interest in illustrator (but that feels like 
cheating).

Jake



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

Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

[Jacob Berv]

This is so cool. Thanks Liam!
J

> On Jun 21, 2018, at 8:22 PM, Liam J. Revell  wrote:
> 
> Dear Theo.
> 
> This can be done with the function rerootingMethod as described here 
> http://blog.phytools.org/2013/04/estimating-ancestral-states-when.html, 
> although it is *very* important to note that rerootingMethod is only valid if 
> the fitted transition model is symmetric (e.g., "ER", "SYM", or any custom 
> model in which i->j == j->i for all i & j). It can also be done with 
> make.simmap for both symmetric & non-symmetric transition models. If you 
> sample enough stochastic maps & then run summary on the set of maps the 
> posterior probabilities will converge on the marginal ancestral states under 
> the defaults for make.simmap, so long as a flat prior is used on internal 
> nodes. This is also on my blog here: 
> http://blog.phytools.org/2013/03/estimating-ancestral-states-when-tips.html.
> 
> 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 6/21/2018 5:32 PM, Théo Léger wrote:
>> Hello,
>> I am working on the phylogeny of a subfamily of moths and I use ace from the 
>> ape-package to reconstruct the ancestral state of a bunch of morphological 
>> characters.
>> I encountered a problem with the few unknown states I have on my matrix 
>> (coded "?", either because material for examination was missing or the state 
>> could not fit in any of the categories): they are treated as other 
>> characters. Is there a way to ignore them? Is there a way to estimate the 
>> state for the species with missing information? I know it is possible to 
>> estimate the state at a tip with missing information in functions like 
>> AncTresh (using bayesian reconstruction) from the phytools package, but I am 
>> not sure if ML reconstruction can do it.
>> Thanks in advance to those who can enlighten me!
>> Cheers,
>> Th�o L�ger
>>  [[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/

___
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] estimate ancestral state with OUwie models

[Jacob Berv]

What about situations in which fossil calibrations are used as priors to inform 
reconstructions for uncalibrated important/interesting nodes? Sure, there is 
great uncertainty, but that doesn’t necessarily imply we should entirely 
abandon hypothesis testing using this approach, does it? I like the idea of 
using priors to constrain different scenarios and then using model testing to 
compare alternative histories (though one’s model ranking may change depending 
on the underlying models too, I guess).

Jake


> On Jun 14, 2018, at 9:32 AM, David Bapst  wrote:
> 
> Simone, Marguerite, others,
> 
> I'll also add that I think there's a great deal to be skeptical of
> ancestral trait reconstruction even when large amounts of fossil data is
> available. You can try the exercise yourself: simulate pure BM on a
> non-ultrametric tree with lots of 'extinct' tips, and you'll still find
> pretty large confidence intervals on the estimates of the trait values.
> What does it mean to do ancestral trait reconstruction, if our calculations
> of uncertainty are that broad?
> 
> That said, in the era of sampled-ancestor phylogenetics for fossil data,
> the ability to examine quantitative support for ancestor-descendant
> relationships among fossil taxa may allow alternative routes to considering
> this issue.
> 
> Cheers,
> -Dave
> 
> On Tue, Jun 12, 2018 at 2:43 AM, Simone Blomberg 
> wrote:
> 
>> I would add an extra caveat to Marguerite’s excellent post: Most
>> researchers work with extant taxa only, ignoring extinction. This causes a
>> massive ascertainment bias, and the character states of the extinct taxa
>> can often be very different to the ancestral state reconstructions,
>> particularly if the evolutionary model is wrong. Eg. there has been an
>> evolutionary trend for example. Ancestral state reconstructions based only
>> on extant taxa should be treated as hypotheses to be tested with fossil
>> data. I wouldn’t rely on them for much more.
>> 
>> Cheers,
>> Simone.
>> 
>> Sent from my iPhone
>> 
>> On 12 Jun 2018, at 4:59 pm, Marguerite Butler 
>> wrote:
>> 
>> Aloha all,
>> 
>> There is no requirement for an ultra metric tree in the formulae reported
>> in Butler-King 2004. Interested investigators should in particular read the
>> supplementary materials where the mathematical details are worked out.
>> 
>> We do generally use ultrametric trees because as comparative biologists,
>> it is more straightforward to think about evolution in units of time rather
>> than in terms of mutational units, etc. However this is by choice, not any
>> methodological limitation.
>> 
>> Once the model parameters are found, the phylogenetic variance-covariance
>> matrix defined by the alpha, thetas, and sigmas can be used to compute
>> ancestral states using a weighted least squares reconstruction method
>> (instead of the typical BM var-cov matrix). The mapping of the alphas,
>> thetas, and sigmas onto the tree are incorporated into this V-COV matrix,
>> so that accounts for the OU model.
>> 
>> NOTES:
>> 1) without knowing why you are doing this, I do feel compelled to warn you
>> that it is unclear why one would want to estimate ancestral states for
>> poorly-fitting models. Be careful!
>> 
>> 2) I hope you realize that ancestral states are in general poorly
>> estimated, even assuming the “correct” model. This is because there is less
>> and less information to anchor the values as you get farther from the tips,
>> similar to the root estimation problem described below.  This issue was
>> clearly exposed in Schluter et al 1997 (and less famously so in Butler and
>> Losos 1997). These depressing results were among the motivations for
>> developing model-fit approaches in the first place.
>> 
>> 3) In 2008/2009 the algorithms in OUCH, SLOUCH, and possibly other methods
>> have changed in the estimation of the value of the root state (X0) which is
>> an internal calculation in fitting the model.  Ho and Ane 2013, and Hansen
>> et al 2008 both reported that the root state X0 is ill-defined (unless
>> there are fossil data to anchor the value). This makes sense intuitively,
>> as all of the information is from the tips, and the root is very far down
>> the tree. A reasonable assumption is that it is distributed according to
>> the stationary distribution of the OU process (X0 ~ N(theta(0),
>> sigma^2/2*alpha) and this assumption is what these methods now employ.
>> 
>> 4) Whatever you end up doing, do check for the robustness of your results
>> with parametric bootstrap on your fitted models (a la Boettinger et al
>> 2012). As many investigators have reported, these parameters can have large
>> confidence intervals, and can covary with one another (being on a
>> likelihood ridge, etc.). But do note that even when parameters may not be
>> uniquely identifiable, it may still be possible to have robust model
>> selection (see Cressler et al 2015).  So perhaps you want to fit ancestral
>> states to see if the different models 

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

[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] Interpretation of standard errors of parameter estimates in OUwie models

[Jacob Berv]

Even if those models may fit the data better, they may not necessarily help 
Rafael determine whether or not the parameter estimates of interest are 
different across regimes (though perhaps BMS might be informative).

Rafael, maybe you could try fixing the ancestral regimes to match most likely 
states for each node from your SIMMAP summary? I wonder if that might decrease 
your ‘uncertainty’ in parameter estimates. 

I don’t have a great answer for your main question though, which is a good one. 

Jake

> On Apr 4, 2018, at 8:59 PM, William Gearty  wrote:
> 
> Hi Rafael,
> 
> Have you tried running the BM1, BMS, or OU1 models? If so, how do the AICc
> values for those compare to those of the OUM model? If not, you should make
> sure you run those.
> If you find that the these models fit your data better, this could indicate
> that there isn't a large different across regimes and an OUM model isn't
> necessarily suitable.
> 
> Best,
> Will
> 
> On Wed, Apr 4, 2018 at 12:30 PM, Rafael S Marcondes > wrote:
> 
>> Dear all,
>> 
>> I'm writing (again!) to ask for help interpreting standard errors of
>> parameter estimates in OUwie models.
>> 
>> I'm using OUwie to examine how the evolution of bird plumage color varies
>> across habitat types (my selective regimes) in a tree of 229 tips. I was
>> hoping to be able to make inferences based on OUMV and OUMVA models, but I
>> was getting nonsensical theta estimates from those. So I've basically given
>> up on them for now.
>> 
>> But even looking at theta estimates from OUM models, I'm getting really
>> large standard errors, often overlapping the estimates from other selective
>> regimes. So I was wondering what that means exactly. How are these erros
>> calculated? How much do high errors it limit the biological inferences I
>> can make? I'm more interested in the relative thetas across regimes than on
>> the exact values (testing the prediction that birds in darker habitats tend
>> to adapt to darker plumages than birds in more illuminated habitats).
>> 
>> I have attached a table averaging parameter estimates and errors from
>> models fitted across a posterior distribution of 100 simmaps for four
>> traits; and one exemplar fitted model from one trait in one of those
>> simmaps.
>> 
>> Thanks a lot for any help,
>> 
>> 
>> *--*
>> *Rafael Sobral Marcondes*
>> PhD Candidate (Systematics, Ecology and Evolution/Ornithology)
>> 
>> Museum of Natural Science 
>> Louisiana State University
>> 119 Foster Hall
>> Baton Rouge, LA 70803, USA
>> 
>> Twitter: @brown_birds 
>> 
>> 
>> ___
>> 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-ph...@r-project.org/
>> 
>> 
> 
> 
> -- 
> William Gearty
> PhD Candidate, Paleobiology
> Department of Geological Sciences
> Stanford School of Earth, Energy & Environmental Sciences
> williamgearty.com
> 
>   [[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/

[R-sig-phylo] plotting geologic time scale on circular tree

[Jacob Berv]

Hi all,
Does anyone know of a function to plot a geologic time scale as a series of 
concentric circles on a circularly plotted tree?

As far as I can tell there are three available functions that can do this on a 
regular cladogram:

axisGeo (phyloch)
geoscale.Phylo  (strap)
geo.legend (phytools)

But none of these works with a circular tree, as far as I can tell. It 
shouldn’t be too hard to code this manually as a series of concentric circles 
(going try that now), but I figured I’d ask here in case someone has already 
done this.

Cheers,
Jake Berv
___
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] Rooting trees in R

[Jacob Berv]

see root() in the R package ape or reroot() in phytools
J

> On May 23, 2017, at 1:22 PM, Valentine Usongo  
> wrote:
> 
> Dear members
> 
> My name is Valentine Usongo and  I am a postdoctoral trainee at the INSPQ 
> Montreal-Quebec. I am working on the use of whole genome sequenced based 
> methods for the detection of foodborne outbreaks in Salmonella. I am still 
> a neophyte in the field of bioinformatics. I just concatenated 14 
> individual trees with the cat function in linux. These trees were in 
> newick format and they were generated using the SNVPhyl pipeline 
> of the National Reference Microbiology Laboratory of Canada. These trees 
> were generated by mapping my Salmonella isolates to a reference genome and 
> the only thing that differs in all the trees is the reference genome since 
> each tree has a different reference and my objective here was to assess 
> the impact of the choice of  reference genome on the ability of the SNV 
> method to separate outbreak from nonoutbreak isolates. My  goal is to 
> compare the topologies of these trees using the Kendall-Colijn (rooted) 
> method. This method requires that trees should be rooted but mine are not 
> rooted. 
> Is it possible to root unrooted trees in R or linux and if so how can one 
> go about it ? The metric requires that the input trees should be in the 
> nexus format and they should be rooted. My trees are in newick format and 
> I converted them to nexus in R but then got stuck because I cannot figure 
> how to root the trees. I would be very grateful for your input. 
> Looking forward to hearing from you. 
> Yours sincerely 
> Valentine Usongo 
>   [[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/

Re: [R-sig-phylo] Graphically comparing multiple trees

[Jacob Berv]

Very interested in this too!
Jake

> On Apr 25, 2017, at 1:20 PM, Vojtěch Zeisek  wrote:
> 
> Hello,
> for comparison of two trees I can use very nice function cophyloplot plotting 
> two trees (left and right) and connecting respective nodes by lines. Very 
> nice 
> and convenient to read. But only for two trees. Displaying multiple trees in 
> multiple comparisons is not very convenient.
> To display dozens to hundreds of trees there is densitree. Also nice, but for 
> this purpose I don't like its display.
> I have several trees (~5) and I wish to compare their topologies, show 
> supports (at least for differing nodes) and highlight differences. I thought 
> about some overlay/parallel plotting (similar to the attached image) where 
> there would be complete topologies displayed and incongruences would be 
> easily 
> visible. It would be probably doable by plotting all separate trees by 
> plot.phylo and then combining and tuning the figure in some vector editor 
> (like Inkscape). But I hope there is some more automated way to do it. :-)
> Sincerely,
> V.
> 
> -- 
> Vojtěch Zeisek
> https://trapa.cz/en/
> 
> Department of Botany, Faculty of Science
> Charles University, Prague, Czech Republic
> https://www.natur.cuni.cz/biology/botany/
> 
> Institute of Botany, Czech Academy of Sciences
> Průhonice, Czech Republic
> http://www.ibot.cas.cz/en/
> ___
> 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/

[R-sig-phylo] phylip file compression ratios

[Jacob Berv]

I noticed today that the compression ratio for an interleaved phylip file (zip 
compressed) was about 84:1, (390MB uncompressed —> 4.6MB compressed) whereas 
the compression ratio for the same data non-interleaved was a much worse 3.4:1 
(390 MB uncompressed —> 113.9 MB). Not knowing much about how zip compression 
actually works - I thought this might be an interesting observation for the 
group…

Best,
Jake Berv
___
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] computing phylogenetic residuals

[Jacob Berv]

Hi all,
I’m trying to compute size corrected phylogenetic residuals on data where the 
independent variable (size) is complete but the dependent variables have 
varying amounts of missing data. phyl.resid() in phytools doesn’t appear to be 
able to handle pairwise dropping (to only consider rows where data is present), 
and outputs only NA when any number of entries in the dependent variables also 
contain NA. Should I just generate new ‘complete’ datasets?

Cheers,
Jake Berv 
___
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] HKY GTR distances

[Jacob Berv]

Can you do this in MEGA?
Jake

> On Feb 1, 2017, at 5:52 PM, kolte...@rub.de  
> wrote:
> 
> Dear Phylothusiasts,
> I need to compare multiple substitution models side-by-side (species 
> clustering stuff by distances only). Unfortunately, I am not aware of an 
> implementation of HKY and GTR distance computations using R. Maybe there is 
> some Github code or something else that I have been missing?
> I do not want to build a phylogenetic tree. I can not use PAUP (>500 big 
> fasta files).
> Any ideas are greatly appreciated.
> Best wishes,
> Andreas
> 
> ___
> 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/

Re: [R-sig-phylo] to log or not to log

[Jacob Berv]

This makes sense to me too… thanks for your thoughts!
Jake

> On Jan 26, 2017, at 7:55 PM, Donald Miles <urosau...@gmail.com> wrote:
> 
> I agree with Ted. The comparison should be how the transformation affects the 
> distribution of the response variable. As Ted mentioned, you should examine 
> the behavior of the residuals to determine whether they conform to the 
> assumptions of the statistical analysis after transformation.
> 
> Donald Miles
> 
> On Thu, Jan 26, 2017 at 7:49 PM, Theodore Garland <theodore.garl...@ucr.edu 
> <mailto:theodore.garl...@ucr.edu>> wrote:
> I don't think you can compare models like this for different transforms of
> the dependent variable.  The likelihood, etc., values are not comparable,
> as your results suggest.  But I am sure someone will correct me if I am
> wrong!
> 
> Assuming this is some sort of regression model (i.e., you have one or more
> independent variables), then you CAN look at residuals to see if they are
> better/worse behaved (e.g., approximately normal, no nasty outliers,
> homoscedasticity).
> 
> Cheers,
> Ted
> 
> 
> On Thu, Jan 26, 2017 at 2:25 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com 
> <mailto:jakeberv.r.sig.ph...@gmail.com>>
> wrote:
> 
> > Dear R-sig-phylo,
> >
> > When analyzing comparative data in a PCM framework, is it be appropriate
> > to use an AIC score to advocate for log transforming input data? ie
> > model(data) vs model(log(data))
> >
> > I’m running a few OUwie models and it’s not entirely clear from the
> > biology whether or not the ‘standard’ log transformation makes sense. In my
> > case, the model(log(data)) results have much lower AIC values relative to
> > the model(data) results.
> >
> > Cheers,
> > Jake Berv
> > ___
> > R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> > <mailto:R-sig-phylo@r-project.org>
> > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> > <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> > Searchable archive at http://www.mail-archive.com/r- 
> > <http://www.mail-archive.com/r->
> > sig-ph...@r-project.org/ <http://sig-ph...@r-project.org/>
> 
> [[alternative HTML version deleted]]
> 
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> <mailto:R-sig-phylo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ 
> <http://www.mail-archive.com/r-sig-phylo@r-project.org/>


[[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] to log or not to log

[Jacob Berv]

Dear R-sig-phylo,

When analyzing comparative data in a PCM framework, is it be appropriate to use 
an AIC score to advocate for log transforming input data? ie model(data) vs 
model(log(data))

I’m running a few OUwie models and it’s not entirely clear from the biology 
whether or not the ‘standard’ log transformation makes sense. In my case, the 
model(log(data)) results have much lower AIC values relative to the model(data) 
results.

Cheers,
Jake Berv
___
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] Comparing support values on different trees

[Jacob Berv]

To clarify - the idea here is that you are asking which clades appear in 
’subordinate’ trees relative to clades that exist in a consensus tree, and then 
interrogating the support values of the shared clades which exist in the 
’subordinate’ trees? So for example, clade A appears in consensus tree X, and 
also appears in gene trees 1-5 - so this will give me the summary of the 
support values for clade A in gene trees 1-5?

Seems like this would be a useful function to interrogate support values among 
clades recapitulated in gene trees relative to a species tree. Would there be 
an analogue for clades that exist in your consensus that don’t exist in 
’subordinate’ trees?

Jake


> On Dec 14, 2016, at 12:41 PM, Keith Barker  wrote:
> 
> Frank:
> 
> You can import all of the trees into one or more multiPhylo objects, then use 
> the ape functions prop.part or prop.clades (depending on what you want to do) 
> to summarize different subsets (e.g., from different analyses). Here is an 
> example with simulated trees:
> 
> x<-rmtree(50,100)
> plot(x[[1]])
> nodelabels(prop.clades(x[[1]],x))
> y<-rep(x[1],50)
> plot(x[[1]])
> nodelabels(prop.clades(x[[1]],c(x,y)))
> 
> The first part just creates a bunch of random trees, so most nodes will only 
> be supported by around 1-4 or so trees. The second part just repeats tree one 
> 50 times, and when you label the nodes with trees x+ tres yy, you get 50 plus 
> the number of trees from part 1. That should give you the idea.
> 
> You can find the shared clades from the "best" trees (if you are doing ML) by 
> first calculating the strict consensus using the consensus function in ape.
> 
> Hope that helps,
> Keith
> 
> On 12/14/16 10:19 AM, Frank T Burbrink wrote:
>> Hello,
>> 
>> I have one question
>> 
>> Is there a method to compare the support values (either bootstraps or Pp) 
>> across all shared clades between two or more different trees having 
>> identical taxa? I believe this method would have to first identify the 
>> shared clades and then determine the measure of support at each shared node.
>> 
>> Thank you,
>> 
>> Frank
>> 
>> Frank T. Burbrink, Ph.D.
>> Associate Curator
>> Department of Herpetology
>> American Museum of Natural History
>> Central Park West at 79th Street
>> New York, NY 10024-5192
>> 
>> fburbr...@amnh.org
>> 
>>  [[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/
> 
> -- 
> F. Keith Barker, Ph.D.
> Associate Professor, Department of Ecology, Evolution and Behavior
> Curator of Genetic Resources, and
> Interim Curator of Birds, Bell Museum of Natural History
> University of Minnesota
> 140 Gortner Laboratory
> 1479 Gortner Ave
> Saint Paul, MN 55108
> 612.624.2737 (phone)
> 612.624.6777 (fax)
> barke...@umn.edu
> http://www.tc.umn.edu/~barke042
> 
> ___
> 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/

Re: [R-sig-phylo] Use of a Phylocom tree for PGLS

[Jacob Berv]

Personally, I don’t think I’d have a problem with this approach (especially 
given the paper Liam cited) given that you are using a phylogeny (a model) to 
test a hypothesis, which is, after all, all we can ever do. You can always do 
more, so any threshold of phylogenetic tree “quality” is going to be somewhat 
arbitrary anyway. I don’t mean to sound defeatist — obviously you should do the 
best job that you can, and it sounds like you have done that given that the 
goal of your research is not to reconstruct a new phylogeny of the clade you’re 
studying. Perhaps you can try to explain this in your rebuttal. 

Or alternatively, use Liam’s awesome code (below) to generate all possible 
phylogenetic trees give your polytomies, and run your PGLS on all of them to 
see if its sensitive to topology. You could probably loop this up somehow 
without much effort. 
http://blog.phytools.org/2016/08/resolving-one-or-more-multifurcations.html?m=1 


Jake



> On Sep 17, 2016, at 6:46 PM, Oscar Valverde  wrote:
> 
> Dear colleages
> 
> I am working on an phylogenetic signal and PGLS analysis using a database
> with values for ~600 plant species. To construct my phylogeny I used the
> backbone Phylomatic supertree (http://phylodiversity.net/phylomatic/) and
> added branch lengths with the based on a fossil calibration for angiosperms
> using the bladj function in phylocom, and resolved polytomies with the
> multi2di function in phylotools.
> 
> Now that I am trying to publish the paper, some reviewers indicated that
> such tree is not suitable for statistical analysis because the level of
> resolution of the tree is too low (to the family level maybe?) and the
> uncertainty is too high to get any reliable result with respect to PGLS or
> phylogenetic signal of the traits. Instead, they suggest I should build my
> own tree based on sequences. Of course this is a major project to undertake
> and in my opinion far from the scope of my study. In fact, this position
> defies the whole reason to have websites like phylomatic where researchers
> can use a reliable resolved phylogenetic tree instead of creating a new one
> every time.
> 
> I would just like to know if the position of my reviewers is a valid one,
> and it the answer is yes, what resource should I use to get a reliable
> phylogenetic tree without making my own version. Thanks in advance for any
> help.
> 
> 
> 
> Oscar Valverde
> Post Doctorate Associate
> International Center on Tropical Botany
> Florida International University
> 
> ‘Anything else you’re interested in is not going to happen if you can’t
> breathe the air and drink the water. Don’t sit this one out. Do something.
> You are by accident of fate alive at an absolutely critical moment in the
> history of our planet.’ ~Carl Sagan
> 
>   [[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 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] suppress axes in phenogram()

[Jacob Berv]

My ugly solution is to comment out lines 195 and 196 in the phenogram code.
Jake

> On Jul 7, 2016, at 5:06 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com> wrote:
> 
> Cool - that works - but now I can’t use axis() to add custom axes…
> Jake
> 
>> On Jul 7, 2016, at 5:03 PM, Liam J. Revell <liam.rev...@umb.edu 
>> <mailto:liam.rev...@umb.edu>> wrote:
>> 
>> par(xaxt="n",yaxt="n")
> 


[[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] suppress axes in phenogram()

[Jacob Berv]

Hi all,
Does anyone know if there is an easy way to suppress the axes in phytools 
phenogram() without altering the code (which is easy to do, but inconvenient 
when using the function in multiple different contexts)? yaxt=’n’ and xaxt=’n’ 
have no effect.

Best,
Jacob Berv
___
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] [function like ratebystate, but for discrete characters?]

[Jacob Berv]

Also see BayesTraits for the bayesian version. 
Jake

> On May 3, 2016, at 8:11 PM, Liam J. Revell  wrote:
> 
> Hi Michael.
> 
> This is essentially the method of Pagel (1994). Though often interpreted as a 
> test for correlated evolution of two binary characters, the model that is 
> actually being fit is one in which the rate of transition 0->1 (or vice 
> versa) in character A is affected by the state of binary character B; or, 
> conversely, if the rate of transition 0->1 (or vice versa) in character B is 
> affected by the state of character A.
> 
> This is implemented in fitPagel of phytools and I believe that this & other 
> more sophisticated flavor are available as part of the package corHMM.
> 
> All the best, Liam
> 
> Liam J. Revell, Associate Professor of Biology
> University of Massachusetts Boston
> web: http://faculty.umb.edu/liam.revell/
> email: liam.rev...@umb.edu
> blog: http://blog.phytools.org
> 
> On 5/3/2016 8:00 PM, Michael Foisy wrote:
>> Hello,
>> 
>> Quick question.  I'm interested in whether the presence/absence of character 
>> 'A' affects the rate of evolution in character 'B'.  Is there a function 
>> that will do this for two or more discrete characters?
>> 
>> I've seen ratebystate {phytools} for continuous data; but my data are all 
>> binary!
>> 
>> 
>> Thank you!
>> , Michael Foisy
>> 
>> MSc. Student
>> University of Toronto
>> michael.fo...@mail.utoronto.ca
>> 
>> 
>>  [[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/

___
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] R-sig-phylo Digest, Vol 95, Issue 6

[Jacob Berv]

There appears to be some issues with ggtree plotting BEAST output - I just 
tried to run through the example and got the following errors: 


source("https://bioconductor.org/biocLite.R;)
biocLite("ggtree")
require("ggtree")

> file <- system.file("extdata/BEAST", "beast_mcc.tree", package="ggtree")
> beast <- read.beast(file)
> beast

Phylogenetic tree with 15 tips and 14 internal nodes.

Tip labels:
A_1995, B_1996, C_1995, D_1987, E_1996, F_1997, ...

Rooted; includes branch lengths.
> str(beast)
List of 23
 $ edge  : int [1:28, 1:2] 16 17 18 18 17 16 19 19 20 20 ...
 $ edge.length   : num [1:28] 3.1 25.44 9.39 6.39 8.82 ...
 $ Nnode : int 14
 $ tip.label : chr [1:15] "A_1995" "B_1996" "C_1995" "D_1987" ...
 $ rate_range_MIN: num [1:14] NA 0.00137 0.00184 0.00119 0.00135 ...
 $ rate_range_MAX: num [1:14] NA 0.00502 0.00422 0.00639 0.00481 ...
 $ length_range_MIN  : num [1:14] NA 0.00603 15.72468 0.42344 0.46641 ...
 $ length_range_MAX  : num [1:14] NA 7.71 36.72 11 3.95 ...
 $ length_95%_HPD_MIN: num [1:14] NA 0.306 22.049 1.88 1.125 ...
 $ length_95%_HPD_MAX: num [1:14] NA 5.18 30.12 6.67 2.93 ...
 $ height: num [1:14] 38.04 34.89 9.47 33.62 31.6 ...
 $ rate  : num [1:14] NA 0.00292 0.00291 0.00291 0.00289 ...
 $ rate_median   : num [1:14] NA 0.00291 0.0029 0.0029 0.00288 ...
 $ length: num [1:14] 0 3.01 25.72 4.35 2.02 ...
 $ length_median : num [1:14] NA 3.01 25.58 4.28 2.01 ...
 $ rate_95%_HPD_MIN  : num [1:14] NA 0.00231 0.00237 0.00228 0.00226 ...
 $ rate_95%_HPD_MAX  : num [1:14] NA 0.00359 0.00348 0.00357 0.0035 ...
 $ height_95%_HPD_MIN: num [1:14] 35.31 31.42 7.32 33.07 30.72 ...
 $ height_95%_HPD_MAX: num [1:14] 40.8 38.7 11.6 34.3 32.6 ...
 $ height_range_MIN  : num [1:14] 33.97 29.2 5.94 33 30.32 ...
 $ height_range_MAX  : num [1:14] 46.9 42.9 15.1 35.4 33.5 ...
 $ height_median : num [1:14] 37.93 34.82 9.39 33.57 31.56 ...
 $ posterior : num [1:14] 1 0.941 1 1 1 ...
 - attr(*, "class")= chr "phylo"
 - attr(*, "origin")= chr 
"/Library/Frameworks/R.framework/Versions/3.2/Resources/library/ggtree/extdata/BEAST/beast_mcc.tree"
> plot(beast, annotation="length_0.95_HPD", branch.length="none") + theme_tree()
NULL
Warning messages:
1: In plot.window(...) : "annotation" is not a graphical parameter
2: In plot.window(...) : "branch.length" is not a graphical parameter
3: In plot.xy(xy, type, ...) : "annotation" is not a graphical parameter
4: In plot.xy(xy, type, ...) :
  "branch.length" is not a graphical parameter
5: In title(...) : "annotation" is not a graphical parameter
6: In title(...) : "branch.length" is not a graphical parameter

So, there is some issue with plotting, and read.beast is also not parsing all 
of the data from internal branches. you can see that it reads data from 14 
branches (I think these are the terminal branches?) but there are actually 28 
total branches in the example. It would be awesome if this actually could 
replace the read.annotated.nexus from OutbreakTools, as this kind of output is 
much easier to workwith.

Jake


> On Dec 10, 2015, at 7:42 AM, Yu, Guangchuang  wrote:
> 
> Dear Roger,
> 
> ggtree can parse BEAST output and visualize BEAST statistics, see the
> document:
> http://www.bioconductor.org/packages/3.2/bioc/vignettes/ggtree/inst/doc/ggtree.html#annotating-tree-with-beast-output
> 
> Bests,
> Guangchuang
> 
> 
> On Thu, Dec 10, 2015 at 7:00 PM,  wrote:
> 
>> Send R-sig-phylo mailing list submissions to
>>r-sig-phylo@r-project.org
>> 
>> To subscribe or unsubscribe via the World Wide Web, visit
>>https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>> or, via email, send a message with subject or body 'help' to
>>r-sig-phylo-requ...@r-project.org
>> 
>> You can reach the person managing the list at
>>r-sig-phylo-ow...@r-project.org
>> 
>> When replying, please edit your Subject line so it is more specific
>> than "Re: Contents of R-sig-phylo digest..."
>> 
>> 
>> Today's Topics:
>> 
>>   1. Re: Plotting sampled-ancestor trees in R (Roger Close)
>>   2. Implementation of Mir et al.'s (2013) tree balanceindex?
>>  (Gabriel Yedid)
>> 
>> 
>> --
>> 
>> Message: 1
>> Date: Wed, 9 Dec 2015 15:50:11 +
>> From: Roger Close 
>> To: David Bapst 
>> Cc: R Sig Phylo Listserv 
>> Subject: Re: [R-sig-phylo] Plotting sampled-ancestor trees in R
>> Message-ID:
>>

Re: [R-sig-phylo] R-sig-phylo Digest, Vol 95, Issue 6

[Jacob Berv]

Whoops! nevermind. It appears read.beast was being masked from read.beast in 
phyloch. Working now.
Jake


> On Dec 11, 2015, at 12:56 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com> 
> wrote:
> 
> There appears to be some issues with ggtree plotting BEAST output - I just 
> tried to run through the example and got the following errors: 
> 
> 
> source("https://bioconductor.org/biocLite.R 
> <https://bioconductor.org/biocLite.R>")
> biocLite("ggtree")
> require("ggtree")
> 
> > file <- system.file("extdata/BEAST", "beast_mcc.tree", package="ggtree")
> > beast <- read.beast(file)
> > beast
> 
> Phylogenetic tree with 15 tips and 14 internal nodes.
> 
> Tip labels:
>   A_1995, B_1996, C_1995, D_1987, E_1996, F_1997, ...
> 
> Rooted; includes branch lengths.
> > str(beast)
> List of 23
>  $ edge  : int [1:28, 1:2] 16 17 18 18 17 16 19 19 20 20 ...
>  $ edge.length   : num [1:28] 3.1 25.44 9.39 6.39 8.82 ...
>  $ Nnode : int 14
>  $ tip.label : chr [1:15] "A_1995" "B_1996" "C_1995" "D_1987" ...
>  $ rate_range_MIN: num [1:14] NA 0.00137 0.00184 0.00119 0.00135 ...
>  $ rate_range_MAX: num [1:14] NA 0.00502 0.00422 0.00639 0.00481 ...
>  $ length_range_MIN  : num [1:14] NA 0.00603 15.72468 0.42344 0.46641 ...
>  $ length_range_MAX  : num [1:14] NA 7.71 36.72 11 3.95 ...
>  $ length_95%_HPD_MIN: num [1:14] NA 0.306 22.049 1.88 1.125 ...
>  $ length_95%_HPD_MAX: num [1:14] NA 5.18 30.12 6.67 2.93 ...
>  $ height: num [1:14] 38.04 34.89 9.47 33.62 31.6 ...
>  $ rate  : num [1:14] NA 0.00292 0.00291 0.00291 0.00289 ...
>  $ rate_median   : num [1:14] NA 0.00291 0.0029 0.0029 0.00288 ...
>  $ length: num [1:14] 0 3.01 25.72 4.35 2.02 ...
>  $ length_median : num [1:14] NA 3.01 25.58 4.28 2.01 ...
>  $ rate_95%_HPD_MIN  : num [1:14] NA 0.00231 0.00237 0.00228 0.00226 ...
>  $ rate_95%_HPD_MAX  : num [1:14] NA 0.00359 0.00348 0.00357 0.0035 ...
>  $ height_95%_HPD_MIN: num [1:14] 35.31 31.42 7.32 33.07 30.72 ...
>  $ height_95%_HPD_MAX: num [1:14] 40.8 38.7 11.6 34.3 32.6 ...
>  $ height_range_MIN  : num [1:14] 33.97 29.2 5.94 33 30.32 ...
>  $ height_range_MAX  : num [1:14] 46.9 42.9 15.1 35.4 33.5 ...
>  $ height_median : num [1:14] 37.93 34.82 9.39 33.57 31.56 ...
>  $ posterior : num [1:14] 1 0.941 1 1 1 ...
>  - attr(*, "class")= chr "phylo"
>  - attr(*, "origin")= chr 
> "/Library/Frameworks/R.framework/Versions/3.2/Resources/library/ggtree/extdata/BEAST/beast_mcc.tree"
> > plot(beast, annotation="length_0.95_HPD", branch.length="none") + 
> > theme_tree()
> NULL
> Warning messages:
> 1: In plot.window(...) : "annotation" is not a graphical parameter
> 2: In plot.window(...) : "branch.length" is not a graphical parameter
> 3: In plot.xy(xy, type, ...) : "annotation" is not a graphical parameter
> 4: In plot.xy(xy, type, ...) :
>   "branch.length" is not a graphical parameter
> 5: In title(...) : "annotation" is not a graphical parameter
> 6: In title(...) : "branch.length" is not a graphical parameter
> 
> So, there is some issue with plotting, and read.beast is also not parsing all 
> of the data from internal branches. you can see that it reads data from 14 
> branches (I think these are the terminal branches?) but there are actually 28 
> total branches in the example. It would be awesome if this actually could 
> replace the read.annotated.nexus from OutbreakTools, as this kind of output 
> is much easier to workwith.
> 
> Jake
> 
> 
>> On Dec 10, 2015, at 7:42 AM, Yu, Guangchuang <g...@connect.hku.hk 
>> <mailto:g...@connect.hku.hk>> wrote:
>> 
>> Dear Roger,
>> 
>> ggtree can parse BEAST output and visualize BEAST statistics, see the
>> document:
>> http://www.bioconductor.org/packages/3.2/bioc/vignettes/ggtree/inst/doc/ggtree.html#annotating-tree-with-beast-output
>>  
>> <http://www.bioconductor.org/packages/3.2/bioc/vignettes/ggtree/inst/doc/ggtree.html#annotating-tree-with-beast-output>
>> 
>> Bests,
>> Guangchuang
>> 
>> 
>> On Thu, Dec 10, 2015 at 7:00 PM, <r-sig-phylo-requ...@r-project.org> wrote:
>> 
>>> Send R-sig-phylo mailing list submissions to
>>>r-sig-phylo@r-project.org
>>> 
>>> To subscribe or unsubscribe via the World Wide Web, visit
>>>https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>>> or, via email, send a message with su

Re: [R-sig-phylo] simulate Ne changes on a given phylogeny

[Jacob Berv]

On a related note- does anyone know of a function or package that can parse and 
store all of the metadata associated with nodes and branches from a BEAST 
output file? read.beast from phyloch will parse all of the data associated with 
internal nodes, but not with terminal branches (for example, substitution rate 
information). I emailed the developer of that tool and he indicated this was by 
design (though he may be able to offer a solution).

Jake

> On Oct 18, 2015, at 12:04 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com> 
> wrote:
> 
> Yes, this has come up in my reading…but there do seem to be situations where 
> ’nearly neutral’ substitutions can be negatively or positively associated 
> with population size:
> 
> From Lanfear 2014 Population size and the rate of evolution:
> 
> "Putting aside variation in the mutation rate, we largely expect the total 
> rate of evolution to be negatively correlated with Ne if slightly deleterious 
> mutations dominate evolution, and to be positively correlated with Ne if 
> advantageous mutations dominate evolution.”
> "On the other hand, any process that leads to an association between Ne and 
> mutation rates will cause a similar association between Ne and neutral and 
> effectively neutral substitution rates. These processes could include effects 
> such as the evolution of mutation rates, and the co-variation of Ne with 
> life-history traits such as generation time”
> 
> I think I need to read more about how to simulate these different kinds of 
> demographic scenarios.
> 
> I suppose the question is, what is the more likely null hypothesis? That 
> mutation rates can change extremely rapidly? Or that demographic fluctuations 
> (pop size, generation time) can induce changes in the detectable substitution 
> rate among lineages. 
> 
> Jake
> 
> 
> 
>> On Oct 17, 2015, at 11:21 PM, Liam J. Revell <liam.rev...@umb.edu> wrote:
>> 
>> Hi Jacob.
>> 
>> Can I add the somewhat boring & probably obvious comment that under the 
>> neutral theory of molecular evolution the substitution rate is independent 
>> of the effective population size?
>> 
>> All the best, Liam
>> 
>> Liam J. Revell, Associate Professor of Biology
>> University of Massachusetts Boston
>> web: http://faculty.umb.edu/liam.revell/
>> email: liam.rev...@umb.edu
>> blog: http://blog.phytools.org
>> 
>> On 10/17/2015 11:01 PM, Jacob Berv wrote:
>>> Hmmm that seems somewhat indirect but might work… I’ll look into that.
>>> 
>>> To give you more information - I’m actually trying to come up with a way to 
>>> test the idea that substitution rate shifts detected with a relaxed 
>>> molecular clock (BEAST) may be driven by changes in effective population 
>>> size. Simulating data for particular scenarios, and then running that 
>>> simulated data through BEAST could be a useful way to test some explicit 
>>> hypotheses I’m interested in. But I have to simulate the data first.
>>> 
>>> Jake
>>> 
>>> 
>>>> On Oct 17, 2015, at 10:40 PM, Brian O'Meara <omeara.br...@gmail.com> wrote:
>>>> 
>>>> Dick Hudson's ms software can simulate gene trees along a species tree or 
>>>> network with migration, changing population size, etc. The package 
>>>> phyclust can call ms. You could then just simulate nucleotides on these 
>>>> gene trees.
>>>> 
>>>> Best,
>>>> Brian
>>>> 
>>>> ___
>>>> Brian O'Meara
>>>> Associate Professor
>>>> Dept. of Ecology & Evolutionary Biology
>>>> U. of Tennessee, Knoxville
>>>> http://www.brianomeara.info <http://www.brianomeara.info/>
>>>> 
>>>> Postdoc collaborators wanted: http://nimbios.org/postdocs/ 
>>>> <http://nimbios.org/postdocs/>
>>>> Calendar: http://www.brianomeara.info/calendars/omeara 
>>>> <http://www.brianomeara.info/calendars/omeara>
>>>> On Sat, Oct 17, 2015 at 10:12 PM, Jacob Berv 
>>>> <jakeberv.r.sig.ph...@gmail.com <mailto:jakeberv.r.sig.ph...@gmail.com>> 
>>>> wrote:
>>>> Dear R-sig-phylo,
>>>> 
>>>> I have a somewhat general simulation question and I was hoping someone on 
>>>> here might have some insight.
>>>> 
>>>> I’m trying to figure out if it’s possible to simulate nucleotide sequence 
>>>> data (an arbitrary number of neutral loci under a multi species coalescent 
>>>> model), on an

Re: [R-sig-phylo] simulate Ne changes on a given phylogeny

[Jacob Berv]

Yes, this has come up in my reading…but there do seem to be situations where 
’nearly neutral’ substitutions can be negatively or positively associated with 
population size:

From Lanfear 2014 Population size and the rate of evolution:

"Putting aside variation in the mutation rate, we largely expect the total rate 
of evolution to be negatively correlated with Ne if slightly deleterious 
mutations dominate evolution, and to be positively correlated with Ne if 
advantageous mutations dominate evolution.”
"On the other hand, any process that leads to an association between Ne and 
mutation rates will cause a similar association between Ne and neutral and 
effectively neutral substitution rates. These processes could include effects 
such as the evolution of mutation rates, and the co-variation of Ne with 
life-history traits such as generation time”

I think I need to read more about how to simulate these different kinds of 
demographic scenarios.

I suppose the question is, what is the more likely null hypothesis? That 
mutation rates can change extremely rapidly? Or that demographic fluctuations 
(pop size, generation time) can induce changes in the detectable substitution 
rate among lineages. 

Jake



> On Oct 17, 2015, at 11:21 PM, Liam J. Revell <liam.rev...@umb.edu> wrote:
> 
> Hi Jacob.
> 
> Can I add the somewhat boring & probably obvious comment that under the 
> neutral theory of molecular evolution the substitution rate is independent of 
> the effective population size?
> 
> All the best, Liam
> 
> Liam J. Revell, Associate Professor of Biology
> University of Massachusetts Boston
> web: http://faculty.umb.edu/liam.revell/
> email: liam.rev...@umb.edu
> blog: http://blog.phytools.org
> 
> On 10/17/2015 11:01 PM, Jacob Berv wrote:
>> Hmmm that seems somewhat indirect but might work… I’ll look into that.
>> 
>> To give you more information - I’m actually trying to come up with a way to 
>> test the idea that substitution rate shifts detected with a relaxed 
>> molecular clock (BEAST) may be driven by changes in effective population 
>> size. Simulating data for particular scenarios, and then running that 
>> simulated data through BEAST could be a useful way to test some explicit 
>> hypotheses I’m interested in. But I have to simulate the data first.
>> 
>> Jake
>> 
>> 
>>> On Oct 17, 2015, at 10:40 PM, Brian O'Meara <omeara.br...@gmail.com> wrote:
>>> 
>>> Dick Hudson's ms software can simulate gene trees along a species tree or 
>>> network with migration, changing population size, etc. The package phyclust 
>>> can call ms. You could then just simulate nucleotides on these gene trees.
>>> 
>>> Best,
>>> Brian
>>> 
>>> ___
>>> Brian O'Meara
>>> Associate Professor
>>> Dept. of Ecology & Evolutionary Biology
>>> U. of Tennessee, Knoxville
>>> http://www.brianomeara.info <http://www.brianomeara.info/>
>>> 
>>> Postdoc collaborators wanted: http://nimbios.org/postdocs/ 
>>> <http://nimbios.org/postdocs/>
>>> Calendar: http://www.brianomeara.info/calendars/omeara 
>>> <http://www.brianomeara.info/calendars/omeara>
>>> On Sat, Oct 17, 2015 at 10:12 PM, Jacob Berv 
>>> <jakeberv.r.sig.ph...@gmail.com <mailto:jakeberv.r.sig.ph...@gmail.com>> 
>>> wrote:
>>> Dear R-sig-phylo,
>>> 
>>> I have a somewhat general simulation question and I was hoping someone on 
>>> here might have some insight.
>>> 
>>> I’m trying to figure out if it’s possible to simulate nucleotide sequence 
>>> data (an arbitrary number of neutral loci under a multi species coalescent 
>>> model), on an ultrametric input topology (where tips represent species), 
>>> with user defined changes in effective population size at the start and end 
>>> of a particular internal branch. In my searching I’ve come across some 
>>> software by Deren Eaton (https://github.com/dereneaton/simLoci 
>>> <https://github.com/dereneaton/simLoci> 
>>> <https://github.com/dereneaton/simLoci 
>>> <https://github.com/dereneaton/simLoci>>) that looks like it might do what 
>>> I want - but I’m not sure. It looks like I can specify migration events 
>>> between taxa, but perhaps not population size changes on internal branches. 
>>> There are many other applications for simulating sequence data but I am not 
>>> familiar with any of them. Any thoughts would be appreciated!
>>> 
>>> Thanks,
>>> 
>>> Jacob Berv
>>> 
>

Re: [R-sig-phylo] simulate Ne changes on a given phylogeny

[Jacob Berv]

Hmmm that seems somewhat indirect but might work… I’ll look into that.

To give you more information - I’m actually trying to come up with a way to 
test the idea that substitution rate shifts detected with a relaxed molecular 
clock (BEAST) may be driven by changes in effective population size. Simulating 
data for particular scenarios, and then running that simulated data through 
BEAST could be a useful way to test some explicit hypotheses I’m interested in. 
But I have to simulate the data first.

Jake


> On Oct 17, 2015, at 10:40 PM, Brian O'Meara <omeara.br...@gmail.com> wrote:
> 
> Dick Hudson's ms software can simulate gene trees along a species tree or 
> network with migration, changing population size, etc. The package phyclust 
> can call ms. You could then just simulate nucleotides on these gene trees. 
> 
> Best,
> Brian
> 
> ___
> Brian O'Meara
> Associate Professor
> Dept. of Ecology & Evolutionary Biology
> U. of Tennessee, Knoxville
> http://www.brianomeara.info <http://www.brianomeara.info/>
> 
> Postdoc collaborators wanted: http://nimbios.org/postdocs/ 
> <http://nimbios.org/postdocs/>
> Calendar: http://www.brianomeara.info/calendars/omeara 
> <http://www.brianomeara.info/calendars/omeara>
> On Sat, Oct 17, 2015 at 10:12 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com 
> <mailto:jakeberv.r.sig.ph...@gmail.com>> wrote:
> Dear R-sig-phylo,
> 
> I have a somewhat general simulation question and I was hoping someone on 
> here might have some insight.
> 
> I’m trying to figure out if it’s possible to simulate nucleotide sequence 
> data (an arbitrary number of neutral loci under a multi species coalescent 
> model), on an ultrametric input topology (where tips represent species), with 
> user defined changes in effective population size at the start and end of a 
> particular internal branch. In my searching I’ve come across some software by 
> Deren Eaton (https://github.com/dereneaton/simLoci 
> <https://github.com/dereneaton/simLoci> 
> <https://github.com/dereneaton/simLoci 
> <https://github.com/dereneaton/simLoci>>) that looks like it might do what I 
> want - but I’m not sure. It looks like I can specify migration events between 
> taxa, but perhaps not population size changes on internal branches. There are 
> many other applications for simulating sequence data but I am not familiar 
> with any of them. Any thoughts would be appreciated!
> 
> Thanks,
> 
> Jacob Berv
> 
> Ph.D. Student
> Lovette Lab
> Cornell Laboratory of Ornithology
> 
> 
> [[alternative HTML version deleted]]
> 
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> <mailto:R-sig-phylo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ 
> <http://www.mail-archive.com/r-sig-phylo@r-project.org/>

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology


[[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] simulate Ne changes on a given phylogeny

[Jacob Berv]

Dear R-sig-phylo,

I have a somewhat general simulation question and I was hoping someone on here 
might have some insight.

I’m trying to figure out if it’s possible to simulate nucleotide sequence data 
(an arbitrary number of neutral loci under a multi species coalescent model), 
on an ultrametric input topology (where tips represent species), with user 
defined changes in effective population size at the start and end of a 
particular internal branch. In my searching I’ve come across some software by 
Deren Eaton (https://github.com/dereneaton/simLoci 
<https://github.com/dereneaton/simLoci>) that looks like it might do what I 
want - but I’m not sure. It looks like I can specify migration events between 
taxa, but perhaps not population size changes on internal branches. There are 
many other applications for simulating sequence data but I am not familiar with 
any of them. Any thoughts would be appreciated!

Thanks,

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology


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

Re: [R-sig-phylo] rotate nodes to match taxonomic order

[Jacob Berv]

The rotateConstr() function does not give me any errors - 

Here is how I was trying to use it:

With a 200 taxon tree (tree A), I rotated ~50 nodes using rotate() to get it 
how I wanted.

With a mostly similar (topologically identical, but different branch lengths, 
so node numbers are not exactly the same as in tree A) 200 taxon tree (tree B), 
I tried:

newtree - rotateConstr(treeB, treeA$tip.label)

this produced a newtree but the ordering looks nothing like treeA - it’s mostly 
jumbled. The root is changed etc… 

Is this the incorrect use of rotateConstr()?

Jake


 On Apr 8, 2015, at 10:57 AM, Emmanuel Paradis emmanuel.para...@ird.fr wrote:
 
 Hi Jacob,
 
 Can you send an example of your output with rotateConstr()? Thanks.
 
 Best,
 
 Emmanuel
 
 Le 08/04/2015 06:41, Jacob Berv a écrit :
 Hi all,
 
 Is there an easy way to get R to automatically rotate the nodes of a 
 phylogeny to match an arbitrary ordering of the tips?
 
 Consider the following two situations -
 
 Situation 1:
 
 Say I have a particular taxonomic order, such as
 
 SpeciesA, SpeciesC, SpeciesB
 
 And I want to rotate the nodes of ((C,B),A) to match it - ie to 
 automatically rotate the nodes to give (A(C,B))
 
 
 Situation 2:
 
 Say I have tree 1 ((C,B),A) and I want to rotate it’s nodes to match the 
 order of tree 2 (A,(B,C))
 
 Currently the only way I know how to accomplish either scenario is to use 
 the ape function rotate() on each relevant node, but this quickly becomes a 
 very tedious task when you have hundreds of nodes to go through and want to 
 achieve a particular ordering.
 
 Any thoughts or tips on how to accomplish either of the two scenarios I 
 describe above in a generalizable way that scales to hundreds of tips/nodes?
 
 There is an example in the rotate() documenatation using rotateConstr():
 # a simple example for rotateConstr:
 A - read.tree(text = ((A,B),(C,D));)
 B - read.tree(text = (((D,C),B),A);)
 B - rotateConstr(B, A$tip.label)
 plot(A); plot(B, d = l)
 
 But this doesn’t seem to work when I try it on my larger trees. Not sure 
 where I’m going wrong…
 
 Best,
 
 Jacob Berv
 
 Ph.D. Student
 Lovette Lab
 Cornell Laboratory of Ornithology
 
 
  [[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/

Re: [R-sig-phylo] multipage pdf of a huge tree

[Jacob Berv]

AFAIK almost any modern browser can render SVG natively - this is not the case 
for PDF, which will often require a plug-in from adobe.
I don’t think I’ve seen this before, but you could probably generate 
web-friendly vector tree figures using SVG - there are options to explore this 
in adobe illustrator if you have it. 

Jake


 On Mar 8, 2015, at 12:06 PM, Jonas Eberle eberle.jo...@gmail.com wrote:
 
 Hi Eike,
 
 I didn't intend to print the tree. I think a pdf is the most convenient
 format for publishing since nearly everybody can easily open it on every
 operating system and it always looks the same. I've never tried to open a
 large tree in an svg in a browser. Maybe that's fine, too.
 
 All the best,
 Jonas
 
   [[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/

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology

___
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] multipage pdf of a huge tree

[Jacob Berv]

Mind=Blown. 

Amazing. Thank you. 

Jake

 On Mar 5, 2015, at 2:55 PM, Liam J. Revell liam.rev...@umb.edu wrote:
 
 Hi Jacob.
 
 This is surprising easy. Check out the solution here: 
 http://blog.phytools.org/2015/03/splitting-tree-across-multiple-pages.html.
 
 You may find that, as in this example, some things like node  edge labels 
 are split across pages - in which case you can fiddle with the split 
 positions  replot the tree until it is the way you want it.
 
 All the best, Liam
 
 Liam J. Revell, Assistant Professor of Biology
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://blog.phytools.org
 
 On 3/5/2015 2:24 PM, Jacob Berv wrote:
 Hi Liam,
 
 Would there be a way to modify this so that annotations like node labels and 
 axes could be applied to both 1/2s of a split tree at the same time? I’ve 
 been struggling with a similar problem using plot.phylo() and we can get our 
 tree to look like what we want on one page, but splitting it up while 
 keeping node and branch annotations etc is problematic.
 
 For example, I am trying to plot a tree and annotate with the following 
 function calls.
 
 #from ape
 plot.phylo(timetree, cex=0.315, label.offset=0.35, edge.width=1.25, 
 edge.color=color)
 
 #from phyloch - adds a geologic time scale
 axisGeo(GTS = strat2012, unit = c(epoch, period), cex = 0.7, col = 
 c(grey80, white), texcol = black, gridty = 2, gridcol = red, 
 ages=TRUE)
 
 #then additional annotations with:
 nodelabels()
 axis()
 
 I want to cut my tree in 1/2 and plot each 1/2 on a separate page, but 
 maintain the axes and annotations. Splitting a tree is great but it’s not 
 super useful unless you can call other functions on the pieces. It would be 
 great to be able to generate axes for each piece of a split tree 
 independently and then annotate each one as if they were a separate plot.
 
 Jake
 
 
 
 On Mar 5, 2015, at 11:46 AM, Jonas Eberle eberle.jo...@gmail.com wrote:
 
 Hi Liam,
 
 perfect! Your right, I didn't realize that some species were duplicated.
 Great to have the option in phytools!
 
 All the best,
 Jonas
 
 
 2015-03-05 17:16 GMT+01:00 Liam J. Revell liam.rev...@umb.edu:
 
 Hi Jonas.
 
 I thought of the same thing. Unfortunately, one thing that I discovered is
 that because R automatically adds 4% or so to the plotting range, you will
 find that species  edges near the bottom of one page will also appear at
 the top of the next page. If you test your code  look for species or edges
 near the bottom or top of each page, you may find what I did.
 
 I have posted a solution to this here: http://blog.phytools.org/2015/
 03/splitting-tree-over-mutiple-plotting.html. It uses plotTree from
 phytools internally, but it could fairly easily be modified to use other
 tree plotters such as plot.phylo or plotSimmap.
 
 All the best, Liam
 
 Liam J. Revell, Assistant Professor of Biology
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://blog.phytools.org
 
 On 3/5/2015 5:19 AM, Jonas Eberle wrote:
 
 MultipagePhyloPlot - function (tree, npages = 5) {
   ntips - length(tree$tip.label)
   tips.per.page - ntips / npages
   y1 - ntips - tips.per.page
   y2 - ntips
   for (i in 1:npages) {
 plot.phylo(tree, y.lim = c(y1, y2))
 y1 - y1 - tips.per.page
 y2 - y2 - tips.per.page
   }
 }
 
 
 
 [[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/
 
 Jacob Berv
 
 Ph.D. Student
 Lovette Lab
 Cornell Laboratory of Ornithology
 

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology

___
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] multipage pdf of a huge tree

[Jacob Berv]

OK - Actually I spoke too soon. That does work for nodelabels() and axis(), but 
it does not work for the axisGeo function in phyloch - this only appears at the 
bottom of the second half of the split plot when I split my tree in 1/2. Is 
there an easy fix to get it to appear on both 1/2s? 


 axisGeo(GTS = strat2012, unit = c(epoch, period), cex = 0.7, col = 
 c(grey80, white), texcol = black, gridty = 2, gridcol = red, 
 ages=TRUE)


Jake


 On Mar 5, 2015, at 2:55 PM, Liam J. Revell liam.rev...@umb.edu wrote:
 
 Hi Jacob.
 
 This is surprising easy. Check out the solution here: 
 http://blog.phytools.org/2015/03/splitting-tree-across-multiple-pages.html.
 
 You may find that, as in this example, some things like node  edge labels 
 are split across pages - in which case you can fiddle with the split 
 positions  replot the tree until it is the way you want it.
 
 All the best, Liam
 
 Liam J. Revell, Assistant Professor of Biology
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://blog.phytools.org
 
 On 3/5/2015 2:24 PM, Jacob Berv wrote:
 Hi Liam,
 
 Would there be a way to modify this so that annotations like node labels and 
 axes could be applied to both 1/2s of a split tree at the same time? I’ve 
 been struggling with a similar problem using plot.phylo() and we can get our 
 tree to look like what we want on one page, but splitting it up while 
 keeping node and branch annotations etc is problematic.
 
 For example, I am trying to plot a tree and annotate with the following 
 function calls.
 
 #from ape
 plot.phylo(timetree, cex=0.315, label.offset=0.35, edge.width=1.25, 
 edge.color=color)
 
 #from phyloch - adds a geologic time scale
 axisGeo(GTS = strat2012, unit = c(epoch, period), cex = 0.7, col = 
 c(grey80, white), texcol = black, gridty = 2, gridcol = red, 
 ages=TRUE)
 
 #then additional annotations with:
 nodelabels()
 axis()
 
 I want to cut my tree in 1/2 and plot each 1/2 on a separate page, but 
 maintain the axes and annotations. Splitting a tree is great but it’s not 
 super useful unless you can call other functions on the pieces. It would be 
 great to be able to generate axes for each piece of a split tree 
 independently and then annotate each one as if they were a separate plot.
 
 Jake
 
 
 
 On Mar 5, 2015, at 11:46 AM, Jonas Eberle eberle.jo...@gmail.com wrote:
 
 Hi Liam,
 
 perfect! Your right, I didn't realize that some species were duplicated.
 Great to have the option in phytools!
 
 All the best,
 Jonas
 
 
 2015-03-05 17:16 GMT+01:00 Liam J. Revell liam.rev...@umb.edu:
 
 Hi Jonas.
 
 I thought of the same thing. Unfortunately, one thing that I discovered is
 that because R automatically adds 4% or so to the plotting range, you will
 find that species  edges near the bottom of one page will also appear at
 the top of the next page. If you test your code  look for species or edges
 near the bottom or top of each page, you may find what I did.
 
 I have posted a solution to this here: http://blog.phytools.org/2015/
 03/splitting-tree-over-mutiple-plotting.html. It uses plotTree from
 phytools internally, but it could fairly easily be modified to use other
 tree plotters such as plot.phylo or plotSimmap.
 
 All the best, Liam
 
 Liam J. Revell, Assistant Professor of Biology
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://blog.phytools.org
 
 On 3/5/2015 5:19 AM, Jonas Eberle wrote:
 
 MultipagePhyloPlot - function (tree, npages = 5) {
   ntips - length(tree$tip.label)
   tips.per.page - ntips / npages
   y1 - ntips - tips.per.page
   y2 - ntips
   for (i in 1:npages) {
 plot.phylo(tree, y.lim = c(y1, y2))
 y1 - y1 - tips.per.page
 y2 - y2 - tips.per.page
   }
 }
 
 
 
 [[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/
 
 Jacob Berv
 
 Ph.D. Student
 Lovette Lab
 Cornell Laboratory of Ornithology
 

___
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] axisChrono and axisPhylo

[Jacob Berv]

Greetings R-sig-phylo,
This is my first post but happy to join this community

I have been having some trouble with the axisChrono of the phyloch package and 
axisPhylo functions of the ape package that I was hoping someone might be able 
to help me with.

When using either function, the spacing between ticks on the scale is too large 
for the figure I’m trying to generate. I have tried passing variations of 
at=c(1,2,3,4,5), or at=seq(0,10,1) etc, but for either function I get the 
following error:

Error in axis(side = side, at = c(maxi - x), labels = abs(x * fact), ...) : 
 formal argument “at matched by multiple actual arguments

Does anyone know how I might be able to specify the tick marks that I want for 
these particular functions?

Best,

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology

___
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] axisChrono and axisPhylo

[Jacob Berv]

Thanks Franz,
I tried modifying the axisChrono function, but when I change the “at=“ argument 
to equal something like “at= c(1,2,3,4,5)” I can get it to list out specific 
tick marks but they are always in increasing order, not reverse order as you 
want with a right facing cladogram. I’m not sure how to reverse the sequence 
within the function.

Jake


Here is the function. 

function (side = 1, unit = Ma, fact = 1, ...) 
{
lastPP - get(last_plot.phylo, envir = .PlotPhyloEnv)
if (lastPP$type %in% c(phylogram, cladogram)) {
if (lastPP$direction %in% c(rightwards, leftwards)) {
x - pretty(lastPP$xx)
if (lastPP$direction == rightwards) 
maxi - max(lastPP$xx)
else {
maxi - min(lastPP$xx)
x - -x
}
}
else {
x - pretty(lastPP$yy)
if (lastPP$direction == upwards) 
maxi - max(lastPP$yy)
else {
maxi - min(lastPP$yy)
x - -x
}
}
}
axis(side = side, at = c(maxi - x), labels = abs(x * fact), 
 ...)
mtext(text = unit, side = side, at = 1.07 * maxi, line = 1, 
  ...)
}





 On Mar 2, 2015, at 5:03 PM, Franz Krah f.k...@mailbox.org wrote:
 
 Hi Jacob,
 
 as far as I can see 'axisPhylo' calculates the at argument of 'axis' 
 itself...
 I think that produces the error if you just put the argument into the 
 'axisPhylo' function. 
 So you have to alter the 'axisPhylo' function so you can define your own 
 'axis' at arguments.
 
 Hope that helps and is correct...
 Best Franz
 
 p.s. that's my first answer here but I feel today is a god day to start 
 (thanks Dave Bapst for encouraging ;-))
 
 Jacob Berv jakeberv.r.sig.ph...@gmail.com hat am 2. März 2015 um 22:38
 geschrieben:
 
 
 Greetings R-sig-phylo,
 This is my first post but happy to join this community
 
 I have been having some trouble with the axisChrono of the phyloch package 
 and
 axisPhylo functions of the ape package that I was hoping someone might be 
 able
 to help me with.
 
 When using either function, the spacing between ticks on the scale is too
 large for the figure I’m trying to generate. I have tried passing variations
 of at=c(1,2,3,4,5), or at=seq(0,10,1) etc, but for either function I get the
 following error:
 
 Error in axis(side = side, at = c(maxi - x), labels = abs(x * fact), ...) : 
 formal argument “at matched by multiple actual arguments
 
 Does anyone know how I might be able to specify the tick marks that I want 
 for
 these particular functions?
 
 Best,
 
 Jacob Berv
 
 Ph.D. Student
 Lovette Lab
 Cornell Laboratory of Ornithology
 
 ___
 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/
 
 MSc. Franz Krah
 Mycology, Phylogenetics, PCM
 http://www.biodiv.wzw.tum.de/index.php?id=18

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology

___
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] axisChrono and axisPhylo

[Jacob Berv]

Excellent! Thanks Santiago. 

In case anyone is interested, I ended up using the following three lines to 
generate an axis with labeled and thick tick marks every 10 and thin 
(unlabeled) ticks every 2. There is probably a more elegant way to do this but 
this works:

axis(3, at=max(branching.times(timetree))-seq(0,70,10), 
labels=seq(0,70,10),cex.axis=0.5, line=-0.95, lwd=2)
axis(3, at=max(branching.times(timetree))-seq(0,70,2),cex.axis=0.25, 
line=-0.95, lwd=0.75, labels=FALSE, tick=TRUE)
mtext(text = Ma, side = 3, at = 1.04 * max(branching.times(timetree)), line = 
-1)


 On Mar 2, 2015, at 11:10 PM, Santiago Claramunt sclaram...@amnh.org wrote:
 
 axis(1, at=max(branching.times(tree))-0:10, labels=0:10)

Jacob Berv

Ph.D. Student
Lovette Lab
Cornell Laboratory of Ornithology


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