Hi - Can you tell me if i'm doing what I think I'm doing....I'm trying to determing if one set of sequences have the same relative branch lengths as another. This falls into the relative ratio test. However I also want to test if allowing for multiple ratios (i.e. some branch are under less constraint than others) is a better fit to the data.

I estimate a synonymous tree with the following substitutions per site:
(((A: 0.002682, B: 0.004542, C: 0.009079): 0.231832, D: 0.134381): 0.127996, E: 0.329363, F: 0.628702);

Now given a sequence I want to compare two sets of constraints:
Tree_13.A.t:=r*.002682;
Tree_13.B.t:=r*.004542;
Tree_13.C.t:=r*.009079;
Tree_13.Node2.t:=r*.231832;
Tree_13.D.t:=r*.134381;
Tree_13.Node1.t:=r*.127996;
Tree_13.E.t:=r*.329393;
Tree_13.F.t:=r*.628702;

Against constraints with sequences A, B, and C under a different relative rate:
Tree_13.A.t:= r1 *.002682;
Tree_13.B.t:= r1 *.004542;
Tree_13.C.t:= r1 *.009079;
Tree_13.Node2.t:=r*.231832;
Tree_13.D.t:=r*.134381;
Tree_13.Node1.t:=r*.127996;
Tree_13.E.t:=r*.329393;
Tree_13.F.t:=r*.628702;

Now I can compare the two with 1 degree of freedom.

Somehow the ratios end up being greater than 1, which is odd because the substitution rates are lower than my synonymous estimates.

I think this is a result of having constraint Tree_13.A.t, which is the number of transitions? rathen than the substitutions/site? Can I constraint the substitution rate?

Anyway, hopefully you can see what I'm trying to do and help me do it correctly.

Thanks for your help.

Scott

Dear Scott,

Firstly, did you, declare r and r1 as global variables?

global r;
global r1;

The code may behave incorrectly if you didn't.

Secondly, the test as set up, will test whether branches A,B,C are evolving at a relative rate different from the rest of the branches, assuming that the relative ratio test hypothesis (using your synonymous tree as a reference) is valid. I think this is what you want to test, so all seems well.

Thirdly, the 'r' (and 'r1') parameter will also include a multiplicative 'branch scaling factor' (since Branch Length = c * t, where c is some constant determined by substitution biases and base composition, and will depend on the model), hence it is not directly comparable to 1. To rescale it after optimization, you can use something like

Code:

rescaled_r = r*Tree_13.A.t/BranchLength (Tree_13,0)/  

which will multiplies r by the reciprocal of 'c' in the BL = c*t relation.

Cheers,
Sergei