[R-sig-phylo] estimating BLs on a fixed topology

2024-03-05 Thread Karla Shikev
Dear all,

I have a large alignment (>3000 species) and a topology and I'd like to
estimate branch lengths on that fixed topology using a GTR+G model, but it
has been challenging computationally using optim.pml in phangorn. Any ideas
about how to speed things up?

Thanks!

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/


Re: [R-sig-phylo] estimating BLs on a fixed topology

2024-03-05 Thread Nathanael Walker-Hale
Hi Karla,

I would recommend using raxml-ng outside of R and then reading the resulting 
tree back in - the command would be something like:

raxml-ng --evaluate --msa [alignment] --tree [tree] --model GTR+G --opt-model 
on --opt-branches on

The tree with the optimised branch lengths will be in 
[alignment].raxml.bestTree.

Cheers,

Nat

From: R-sig-phylo  on behalf of Karla Shikev 

Sent: Tuesday, March 5, 2024 7:17:00 PM
To: R Sig Phylo Listserv 
Subject: [R-sig-phylo] estimating BLs on a fixed topology

Dear all,

I have a large alignment (>3000 species) and a topology and I'd like to
estimate branch lengths on that fixed topology using a GTR+G model, but it
has been challenging computationally using optim.pml in phangorn. Any ideas
about how to speed things up?

Thanks!

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] chronos ape substitution rate strict clock

2024-03-05 Thread Vincenzo Ellis
Dear Klaus,

Thank you for your response. I was wondering if you could expand on how to
do this with phangorn::pml_bb(). Here's a small tree I made in Raxml with
the cytochrome b gene of a few species of birds:
library(ape)
tree <- read.tree(text =
"(Nesillas_typica:0.04793543089316195,(((A_schoenobaenus:0.043980600360374364,A_melanopogon:0.06481521246716708):0.022244127250028972,(((I_caligata:0.08810019775981455,I_pallida:0.062476654152826855):0.017350459234948976,(((A_newtoni:0.033344055105229795,A_sechellensis:0.01792648331103279):0.021283417245643144,(A_rufescens:0.05238658747602984,A_gracilirostris:0.044877312044462714):0.013125661916305195):0.022403450749406223,((A_arundinaceus:0.04093765705229456,A_stentoreus:0.01826146619561772):0.015107695382970854,A_griseldis:0.0972067122912957):0.011211088579910683):0.02819743631258305):0.007869067285886489,(H_icterina:0.058362559770475525,H_polyglotta:0.033259000190850735):0.05342630181117944):0.005148698839842897):0.006239321505491391,((A_palustris:0.05185723350471064,(A_scirpaceus:0.014053067606248151,A_baeticatus:0.008365145972302118):0.045911856091680914):0.023773812432215297,(A_dumetorum:0.08560190754497965,A_agricola:0.07765117111045171):0.01):0.013260470881828468):0.04793543089316195);")

If we assume that the divergence rate for cytochrome b for these species
is 2.21% per million years, how could we rescale the tree with a strict
clock? Do we need to start with a DNA alignment and re-estimate the tree
topology or can this be done just with the ML tree I already estimated and
the rate?

Thank you,

Vincenzo

On Mon, Mar 4, 2024 at 7:57 PM Klaus Schliep 
wrote:

> Dear Vincenzo,
> For ML estimates the edge length are the expected number of substitutions
> per site, which depends on the product of rate and time. So you need either
> a rate estimate or calibration points to estimate the time. If you divide
> the edge length by the rate the edges should be proportional to time.
> With pml_bb you can estimate ultrametric or tip dated trees with a strict
> molecular clock model.
> Kind regards,
> Klaus
>
>
>
>
> Von meinem iPhone gesendet
>
> > Am 04.03.2024 um 21:06 schrieb Vincenzo Ellis :
> >
> > Dear Liam,
> >
> > Thank you for your answer. Yes, I have a hypothesized average clock rate
> > and no explicit calibration points.
> >
> > If I use ape::chronos() with default values it gives me a depth of 1 at
> the
> > root. So I suppose to rescale the branch lengths I just need to multiply
> > all of the branch lengths by a value (i.e., tree$edge.length * value), is
> > that right?
> >
> > Thanks again,
> >
> > Vincenzo
> >
> >> On Mon, Mar 4, 2024 at 11:43 AM Liam J. Revell 
> wrote:
> >>
> >> Dear Vincenzo.
> >>
> >> If I understand your problem problem, you do not have any explicit
> >> calibration points -- but you have a hypothesized average clock rate?
> >>
> >> If so, then you can obtain an ultrametric tree from* ape::chronos* for
> >> any value of the smoothing parameter (*lambda*) and then simply re-scale
> >> it to have the desired total depth (based on your hypothesized clock
> rate).
> >>
> >> To choose a "correct" value of *lambda* one can use cross-validation as
> >> described in Sanderson (2002;
> doi:10.1093/oxfordjournals.molbev.a003974).
> >>
> >> Others should feel welcome to weigh in if this is not right.
> >>
> >> All the best, Liam
> >> Liam J. Revell
> >> Professor of Biology, University of Massachusetts Boston
> >> Web: http://faculty.umb.edu/liam.revell/
> >> Book: Phylogenetic Comparative Methods in R
> >> <
> https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r>
> >> (*Princeton University Press*, 2022)
> >>
> >>
> >> On 3/4/2024 11:09 AM, Vincenzo Ellis wrote:
> >>
> >> [You don't often get email from vael...@udel.edu. Learn why this is
> important at https://aka.ms/LearnAboutSenderIdentification ]
> >>
> >> CAUTION: EXTERNAL SENDER
> >>
> >> Dear Emmanuel,
> >>
> >> Thank you very much for your response. I cannot see how to provide the
> >> substitution rate to the phangorn::pml_bb() function, but I was looking
> at
> >> the ape::node.dating() function and it appears that I could provide the
> >> substitution rate to the "mu" argument and then set the "node.dates"
> >> argument to NA or zero for all tips (I'm not sure if NA or zero would be
> >> preferable to force the tips to all be from a single time point). Do you
> >> think that would work? I'm not sure how to make ape::node.dating()
> accept a
> >> substitution rate rather than try to estimate one. Maybe an option
> could be
> >> added to allow mu to equal a user-specified number rather than the
> output
> >> of ape::estimate.mu()?
> >>
> >> Another option might be to calculate an estimated age for every node
> >> connecting sister taxa in the tree by converting the genetic distances
> >> between sister pairs to divergence times using the substitution rate and
> >> then use those as priors in ape::chronos(). I suppose I could also apply
> >> that logic t