Hi - I'm having trouble getting started with an analysis I'd like to do. I've searched through the discussion group and haven't really found an answer, so hopefully you can help me.

I want to compare the likelihood that a noncoding sequence is evolving under a neutral model with the likelhood that the same sequence is constrained and evolving slower than neutral.

I've estimated the neutral model (using concatenated 4-fold degenerate sites) using HYPHY and get a tree that looks like:

((A:0.22, B:0.12)Node1:0.12,C:0.30,D:0.54);
TVTS = .209335
So using an HKY model, I'd like to fix the branch lengths and TVTS ratio, and calculate the likelihood of another sequence file using those fixed parameters.

Then I'd like to test if allowing for a contrained branch length paramenter (that is the keeping the branch lengths proportional, but scalling them all by some fraction) provides a better fit to the data.

This seems like a simple analysis, but I'm having a hard time getting it to work.

Do you have example code that might help?

Thanks,
Scott

ok, so using old posts to this discussion group I think I figured this out. Sorry to make you waste your time.

My problem was that I was building a neutral tree, but then optimizing the likelihood function, so it changed the branch lengths (duh.) So if I fix the branch lengths, the unoptimized likelihood is what I want.

I then build a second tree, constraining the branch lengths by a variable R. And get the likelihood of this constrained model, and output the relative constraint.

Here's what I have, can you confirm that this is doing what I think it's doing?

global TvTs:= 0.205; /* the global transition/transversion ratio */

HKYM =  [*,TvTs*t,t,TvTs*t}
   {TvTs*t,*,TvTs*t,t}
   {t,TvTs*t,*,TvTs*t}
   {TvTs*t,t,TvTs*t,*];


eqf = [0.32}{0.18}{0.18}{0.32];


Model  HKY = (HKYM,eqf,1);

ACCEPT_ROOTED_TREES = 1;

t    = 1;  

/* it is important to set the branch length parameter to 1, so that the rate matrix is evaluated properly for the function,
otherwise the scaling factor will be of by a factor of 1/t */

scalingB = computeScalingFactorB (HKYM, eqf);

Tree L_1 = (((Scer:0.23,Spar:0.15):0.121206,Smik:0.32):0.12,Skud:0.329292,Sbay:0.44);
UseModel(HKY);
Tree L_2 = (((Scer:0.23,Spar:0.15):0.121206,Smik:0.32):0.12,Skud:0.329292,Sbay:0.44);


fprintf (stdout, "\nBranch scaling factor computed to be ", scalingB,"\n\n");

/* get all branch names and branch lengths*/

branchNames  = BranchName  (L_1,-1);
branchLengths = BranchLength (L_1,-1);
global R = 1;
for (k = 0; k < Columns (branchNames)-1; k=k+1)
{
ExecuteCommands ("L_1."+branchNames[k]+".t=L_1."+branchNames[k]+".t/scalingB;");
}
branchLengths2 = BranchLength (L_1,-1);
for(k=0; k < Columns(branchNames)-1; k=k+1)
{
ExecuteCommands ("L_2."+branchNames[k]+".t := R * "+branchLengths2[k]+";");

}

characters = ["A","C","G","T"}
   {"1","","",""];

DataSet B = ReadDataFile("temp.fa");
DataSetFilter partition1 = CreateFilter(B,1);

LikelihoodFunction neutral_LF = (partition1,L_1);
LikelihoodFunction constrained_LF = (partition1,L_2);
     
fprintf(stdout, "\nunoptimized\n", neutral_LF,"\n");
fprintf(stdout, "\nunoptimized constrained\n", constrained_LF,"\n");

Optimize(res,neutral_LF);
Optimize(res2,constrained_LF);

fprintf(stdout, "\nneutral optimized\n",neutral_LF);
fprintf(stdout, "\nconstrained opt.\n",constrained_LF);
rescaled_R = R/(L_2.Scer.t/BranchLength(L_2,0));
fprintf(stdout, "\nrescaled R = ", rescaled_R);

Dear Scott,

Looks OK to me. Let me know if you run into any problems.

Cheers,
Sergei