hi,

You guys have built some very nice tools, thanks.

I have a particular problem I'm writing a batch file for, and I'd like to check if I'm doing what I intend to be doing.

I have two corresponding datasets from noncoding sequence. One is putatively neutral, and the other putatively functional. (for example they could be an ancestral repeat, and a putative regulatory element which are nearby each other). I want to test if the substitution rate for a particular branch (branch A) in the functional is higher than expected relative to the neutral. I define "expected" based on information external to this data--say the ratio of subsitution rates genome wide.

For a 3 species tree ((A,B)C), I put each type of sequence in its own partition, and set up a joint likelyhood function.

LikelihoodFunction  theLnLik = (neutFilt, neutTree, consFilt, consTree);

Then to do a likelyhood ratio test, I get the likelyhood under a constrained model. (I'm using HKY85 by the way).

global ss=0;
consTree.A.trst:=ss*0.26*neutTree.A.trst;
consTree.A.trvs:=ss*0.24*neutTree.A.trvs;
consTree.B.trst:=ss*0.68*neutTree.B.trst;
consTree.B.trvs:=ss*0.69*neutTree.B.trvs;
consTree.C.trst:=ss*0.36*neutTree.C.trst;
consTree.C.trvs:=ss*0.35*neutTree.C.trvs;

In the second statement, the number 0.26 is a branch specific scaling factor (obtained externally to this data) giving the average ratio between neutral and conserved for branch A transitions. ss is a scaling factor for the level of constraint of this particular element in all branches of the tree.

I then compare the likelihood for this model with that for a model with one additional parameter. This parameter allows branch A to go more quickly.

global ss=0;
global RR=1;
RR>1;
consTree.A.trst:=ss*RR*0.26*neutTree.A.trst;
consTree.A.trvs:=ss*RR*0.24*neutTree.A.trvs;
consTree.B.trst:=ss*0.68*neutTree.B.trst;
consTree.B.trvs:=ss*0.69*neutTree.B.trvs;
consTree.C.trst:=ss*0.36*neutTree.C.trst;
consTree.C.trvs:=ss*0.35*neutTree.C.trvs;

Is there some problem with my setting RR>1 like this? Are there other problems with how I have written it?

What I want to be testing is whether there has been a relaxation of contraint in branch A.

thanks much,
Eliot

Dear Eliot,

Everything looks OK except that the one sided constraint on 'R' should read:

Code (]
global R:>1;
[/code):

Also, you you should then test the hypothesis that R!=1 by using a mixture 0.5(chi^2_1+chi^2_0), i.e. you if you have the likelihood ratio statistic LR=(2*(LogL_alt - LogL_null)), the p-value will be [code] pValue = 0.5(1-CChi2(LR,1));

Cheers,
Sergei

thanks a lot for the help!