Branch length scaling

[Ugo Bastolla]

Hello Alexandros, Olexsey and all staff,

first of all thank you very much for developing and maintaining raxml, great job. I have a question about branch length scaling. I am performing partitioned analysis with substitution matrices computed for each protein site by my program SSCPE, which takes into account structural constraints on evolution. Since raxml normalizes the substitution rates internally, Olexsey suggested to me some months ago to input a branch length parameter to each partition with the instruction +BU{parameter}.

I try as BU parameter both the site-specific rate and its inverse. I then run raxml with the option --brlen scaled --opt-model off. How is this instruction actually executed by raxml? I asked it long ago, but perhaps not through the proper channel, and I am realizing that I am not completely sure. Reading the manual, I interpreted that raxml computes a unique set of branch lengths for all partitions and then rescales the branch lengths of each partition and recomputes the likelihood. Both the final likelihood and the final tree are not too much dependent on whether I use rate or 1/rate as the BU parameter, and they change little if I rescale all partitions with a constant parameter (say +BU{2}), therefore I think that raxml must do some internal rescaling for avoiding that the likelihood changes much with a constant scale, but I am not very sure on which computation it performs, can you please help me?

Thank you very much,

Ugo

[Oleksiy Kozlov]

Hi Ugo,

> Hello Alexandros, Olexsey and all staff,
> first of all thank you very much for developing and maintaining raxml, great job.

thank you!

> I try as BU parameter both the site-specific rate and its inverse. I then run raxml with the option
> --brlen scaled --opt-model off. How is this instruction actually executed by raxml? I asked it long
> ago, but perhaps not through the proper channel, and I am realizing that I am not completely sure.
> Reading the manual, I interpreted that raxml computes a unique set of branch lengths for all
> partitions and then rescales the branch lengths of each partition and recomputes the likelihood.

not quite, "scaled" means that one global set of branch lengths is computed, and then multiplied
with individual scaling factors for every partition. Please see here:

https://github.com/amkozlov/raxml-ng/wiki/Tutorial#branch-length-linkage

> Both the final likelihood and the final tree are not too much dependent on whether I use rate or
> 1/rate as the BU parameter, and they change little if I rescale all partitions with a constant
> parameter (say +BU{2}), therefore I think that raxml must do some internal rescaling for avoiding
> that the likelihood changes much with a constant scale, but I am not very sure on which computation

Yes, branch length scalers are normalized such as the weighted (w.r.t. number of sites) average of
scalers is 1.0. You can see the respective code here:

https://codeberg.org/togkousa/coraxlib/src/branch/dev/src/corax/tree/treeinfo.c#L1154

This normalization will happen during branch length optimization, so technically you can disable it
with "--opt-branches off".

Hope this helps,
Oleksiy

[Ugo Bastolla]

Thank you very much for the fast answer Oleksiy!
My problems is  how does raxml rescale the branch lengths of each
partition? Does it multiply or divide times the parameter +BU{} that is
input for each partition?
Another question is whether it is possible to obtain in output more
mLTrees than just one per each starting tree, which would be useful for
post-processing
Best,

Ugo

[Oleksiy Kozlov]

> Thank you very much for the fast answer Oleksiy!
> My problems is  how does raxml rescale the branch lengths of each partition? Does it multiply or
> divide times the parameter +BU{} that is input for each partition?

global branch length is multiplied with the partition-specific branch length scaler, see:

https://codeberg.org/togkousa/coraxlib/src/branch/dev/src/corax/tree/treeinfo.c#L1127

> Another question is whether it is possible to obtain in output more mLTrees than just one per each
> starting tree, which would be useful for post-processing

Not directly, you have two option:
- modify source code
- track&copy $PREFIX.raxml.lastTree.TMP which contains the best tree found so far

Best,
Oleksiy