Hello,

I'm trying to use the branch lengths specified in the tree file for the estimation of site specific substitution rates while estimating the model parameters (actually, a GTR) from the data. How do I have to proceed to constrain the branch lengths before the likelihood optimization?
Thanks for your help!

Seraina

Hi Wayne,
thanks for your reply. Actually, I've already used (and modified for my purpose) the batch file SiteRates2.bf. However, when it comes to the step "Optimize (resNull,lf)", not only the model parameters, but also the branch lengths of my input tree are modified. So what I'm looking for is a way to constrain the optimization procedure to only estimate the parameters of the GTR model from the data, but not modifying the branch lengths. Is this possible?
Thanks for the help!
Seraina

Hi Seraina, I'll try answer this post and the one to the feedback board here. so yes there is a difference between SiteRates.bf and SiteRates2.bf.

SiteRates.bf will first fit a model of the entire alignment in which it estimates both branch length parameters and model parameters (AC, GC etc depending on your model). In the site wise model, the model parameters are constrained to those estimated across the alignment, and the branch lengths (per site) are set to someFactorX*alignmentWideBranchLength. So the site wise analysis returns a substitution rate at a site relative to the average across all sites. It does not fit independent model parameters for each site (AC, GT, etc are all shared across all sites).

SiteRates2.bf will do the same, but the BranchLengths in someFactorX*alignmentWideBranchLength are taken from file. This would be useful if you had some independent estimates of length of branches in units of time.

to answer your direct ?

ReplicateConstraint ("this1.?.?:=this2.?.?__",givenTree,givenTree)

in the SiteRates2.bf file will set the givenTree to exact values from the tree read from file. These will apply to all parameters associated with branches, but not with model parameters that are global across all branches.

There was a problem with the previous SiteRates2.bf in that not all the someFactorX's were being applied. I've attached an edit that will also go into the next release. If what I've described above is what you want  then use this code. But from your later posts I think you actually want independent model parameters (GTR) for each site.

hope this helps somewhat ./w

Dear Seraina,

Back in the day, I had a very long discussion with someone else who wanted to do a similar thing: Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login

There are two salient points to be carried over to your case

1). Branch lengths can mean two things: a priori measures (i.e. millions of years since divergence) or [what is usually meant for model fitting] genetic distance. The latter is a product of evolutionary time and evolutionary rate. Because for a single gene, the evolutionary time is presumably the same for every site, if you hold branch lengths (and here I take them to mean genetic distances) the same for every site, that means that the rate is the same for every site. 

In your input tree, what do the branch lengths mean?

2). The reason HyPhy gives you different estimates for 20 sites vs the total alignment is because it estimates model parameters (e.g. AT rate and the base frequencies) either from the first 20 sites, or the whole alignment -- giving different values. Those values will then influence the estimates of site-specific rates. The only way to remove this effect, is to specify the substitution model (base frequencies and substitution rates a priori). I am quite certain this is what DNARates used to do (with F84 model supplied by the user and base frequencies are set to 0.25 or prompted for).

SiteRates2.bf can be modified for this purpose (take a look at the files from the thread I quoted). Let me know if this helps.

Sergei

Dear Seraina,

Code:

x:=x__;
 

can be read as 'constrain x to its current value'. The '__' postfix is the evaluator operation that replaces the variable (x) with its current value before the formula is evaluated.

Cheers,
Sergei