Hi,

I am working with site-specific models and then I have to build a likelihood function containing the n positions in an alignment like this:

LikelihoodFunction lk = (class0,tree0, class1,tree1,class2, tree2....classn,treen);

I have tried using ExecuteCommands as:

ExecuteCommands("LikelihoodFunction lk = (class0,tree0,");
      for(i=1;i<29;i=i+1){
       ExecuteCommands("class"+i+",tree"+i+", ");
      }
       ExecuteCommands("class29,tree29); ");

for an alignment with n=30, but it doesn't work.
Could you please help me with this script?
Thank you very much.

gustavo

Dear Sergei,
Thank you very much for your help!

I have still some problems using site-specific substitution matrices. In the script that follows I would like to use a gamma distribution with the use of site-specific substitution matrices (however here in the example I use JTT model):

    lk=0;
    global alpha = .5;
    alpha:>0.01;alpha:<100;
    category c = (4, EQUAL, MEAN, GammaDist(_x_,alpha,alpha),CGammaDist(_x_,alpha,alpha),0,1e25,CGammaDist(_x_,alp
ha+1,alpha));

    for(i=0;i<len;i=i+1){
       ExecuteCommands("DataSetFilter class"+i+" = CreateFilter (data,1,\""+i+"\",\"\");");
       obsFreqs={20,1};
       obsFreqs=0;
       ExecuteCommands("HarvestFrequencies (obsFreqs"+i+", class"+i+", 1, 1, 1);");
       ExecuteCommands("#include Jones2.mdl\";");
       ExecuteCommands("Model JTT"+i+" = (jones, obsFreqs"+i+",1);");
       ExecuteCommands("Tree tree"+i+" = treestring;");
    }

lfDefCommand = "";
lfDefCommand * "LikelihoodFunction lk = (class0,tree0,";
for(i=1;i<len-1;i=i+1){
        lfDefCommand * ("class"+i+",tree"+i+", ");
}
i=len-1;
lfDefCommand * ("class"+i+",tree"+i+"); ");
lfDefCommand * 0;
fprintf (stdout, "\nExecuting: ", lfDefCommand, "\n");
ExecuteCommands (lfDefCommand);

    for(i=1;i<len;i=i+1){
       ExecuteCommands("ReplicateConstraint (\"this1.?.a:=this2.?.a\",tree"+i+",tree0);");
    }

    Optimize(params,lk);

    fprintf(prueba1.out,"\n jones Log(L) + F + gamma con constraints lk = ", params[1][0]);

The problem with the script is that never ends! Would you please check it?

Additionally, I would like to ask you why is so important to optimize the time as a common variable for all the positions (using ReplicateConstraint in the script). In this way it is true that you have a given branch length common for all the positions for a given topology. However, if you optimize the time for each position and then derives the total likelihood  adding the log lik, the value you obtained is much better than using the time as a common variable among positions.

thank you very much again and merry christmas!

gustavo

Dear Sergei,
After 6 hs of calculation using the script:

REPLACE_TREE_STRUCTURE = 1;

DataSet data = ReadDataFile ("seqfile");

file1 = "treefile";
fscanf (file1,"Tree",tree);
treestring=""+tree;

file2 = "len";
fscanf (file2,"Number",len);

    lk=0;
    global alpha = .5;
    alpha:>0.01;alpha:<100;
    category c = (4, EQUAL, MEAN, GammaDist(_x_,alpha,alpha),CGammaDist(_x_,alpha,alpha),0,1e25,CGammaDist(_x_,alp
ha+1,alpha));

    for(i=0;i<len;i=i+1){
       ExecuteCommands("DataSetFilter class"+i+" = CreateFilter (data,1,\""+i+"\",\"\");");
       obsFreqs={20,1};
       obsFreqs=0;
       ExecuteCommands("HarvestFrequencies (obsFreqs"+i+", class"+i+", 1, 1, 1);");
       ExecuteCommands("#include \"Jones2.mdl\";");
       ExecuteCommands("Model JTT"+i+" = (jones, obsFreqs"+i+",1);");
       ExecuteCommands("Tree tree"+i+" = treestring;");
    }

lfDefCommand = "";
lfDefCommand * "LikelihoodFunction lk = (class0,tree0,";
for(i=1;i<len-1;i=i+1){
        lfDefCommand * ("class"+i+",tree"+i+", ");
}
i=len-1;
lfDefCommand * ("class"+i+",tree"+i+"); ");
lfDefCommand * 0;
fprintf (stdout, "\nExecuting: ", lfDefCommand, "\n");
ExecuteCommands (lfDefCommand);

    for(i=1;i<len;i=i+1){
       ExecuteCommands("ReplicateConstraint (\"this1.?.a:=this2.?.a\",tree"+i+",tree0);");
    }

LIKELIHOOD_FUNCTION_OUTPUT=7;
fprintf ("forsergei.out",CLEAR_FILE,lk);

    Optimize(params,lk);

    fprintf(prueba1.out,"\n jones Log(L) + F + gamma con constraints lk = ", params[1][0]);


the file "forsergei.out" contains a "0"!

Thanks very much for the answer about the estimation site-by-site.
best regards,

gustavo

hi sergei,
here are the files.
regards,

gustavo