Hi Sergei, I am trying to code a generic branch length scaling function that will work with all rate matrices, including non-reversible matrices. The code is pasted below:


Code (]
function BranchScaling ( rateMatrix, freqVector, reverse, useFreqVector )
{
     mSize = Rows ( freqVector );
     stringB = "";
     stringB * 128;
     if ( reverse ) {
           for ( i = 0; i < mSize-1; i = i + 1 ) {
                 for ( j = i + 1; j < mSize; j = j + 1 ) {
                       if ( useFreqVector ) {
                             if ( ( i == 0 ) && ( j == 1 ) ) {
                                   stringB = "2*(" + rateMatrix [i):

[j] + "*" + freqVector[i] + "*" + freqVector[j];                              }                              else {                                    stringB = stringB + "+" + rateMatrix[i][j] + "*" + freqVector[i] + "*" + freqVector[j];                              }                        }                        else {                              if ( ( i == 0 ) && ( j == 1 ) ) {                                    stringB = "2*(" + rateMatrix[i][j];                              }                              else {                                    stringB = stringB + "+" + rateMatrix[i][j];                              }                        }                  }                  }      }      else {            for ( i = 0; i < mSize-1; i = i + 1 ) {                  for ( j = i + 1; j < mSize; j = j + 1 ) {                        if ( useFreqVector ) {                              if ( ( i == 0 ) && ( j == 1 ) ) {                                    stringB = "(" + rateMatrix[i][j] + "*" + freqVector[i] + "*" + freqVector[j];                                    stringB = stringB + "+" + rateMatrix[j][i] + "*" + freqVector[j]+ "*" + freqVector[i];                              }                              else {                                    stringB = stringB + "+" + rateMatrix[i][j] + "*" + freqVector[i] + "*" + freqVector[j];                                    stringB = stringB + "+" + rateMatrix[j][i] + "*" + freqVector[j] + "*" + freqVector[i];                              }                        }                        else {                              if ( ( i == 0 ) && ( j == 1 ) ) {                                    stringB = "(" + rateMatrix[i][j];                                    stringB = stringB + "+" + rateMatrix[j][i];                              }                              else {                                    stringB = stringB + "+" + rateMatrix[i][j];                                    stringB = stringB + "+" + rateMatrix[j][i];                              }                        }                  }            }      }      stringB = stringB + ");";      /* fprintf ( stdout, stringB, "\n" ); */      stringB * 0;      return stringB; }

The idea would be to call the function with something like:

global B := BranchScaling ( NucRateMatrix, vectorOfFreqs, 1, 0);

This is currently not working, since rateMatrix[i][j] currenlty passes a value to the stringB and not the actual components.

I get:

global B := 2*(0.143358*0.398282*0.171019+...

instead of:

global B := 2*(kappa*0.398282*0.171019+...

As I understand it B should be optimised along with kappa. I have tried various options with no success as yet. It should be possible to do since a LIKELIHOOD_FUNCTION_OUTPUT value of 6 returns the components of rate matrices:

eg:
NucRateMatrix[0][1]:=kappa*t;
NucRateMatrix[0][2]:=t;
NucRateMatrix[0][3]:=kappa*t;
NucRateMatrix[1][0]:=kappa*t;
NucRateMatrix[1][2]:=kappa*t;

Any help/suggestions would be appreciated. I think this would be a neat function to have - especially when scaling more complicated codon models to nucleotide branch lengths. thanks ./w

Oops, the "stringB * 0;" in the third line of code shouldn't be there. ./w

Dear Wayne,

I have been meaning to build in a function for this, since the usage is quite common and many of the expressions can be greatly simplified internally. Let me put that into the next release...

Cheers,
Sergei

Dear Wayne,

As of today's build calling

[code]
GetString(stringID,model_id,-1)
[/code]

will store the expression for the branch length of model model_id into string variable stringID.

Cheers,
Sergei

thanks ./w

Hi Sergei, a quick question regarding GetString. The string returned includes t which as i understand shouldn't be included in the scaling factor. Any suggestions on how to remove t? some test code (modified from previous posts) is attached. cheers ./w


Code (]
ACCEPT_BRANCH_LENGTHS = 1;
ACCEPT_ROOTED_TREES = 1;

DataSet ds = ReadDataFile ( "10taxa.nex" );
fscanf ( "10taxa_phy.tre", "String", TREE );
DataSetFilter filteredData = CreateFilter ( ds, 1 );

UseModel (USE_NO_MODEL);  
Tree dummy = TREE;
branchNames     = BranchName   (dummy,-1);
branchLengthsNow  = BranchLength (dummy,-1);

global AC = 0.1;
global AT = 0.1;
global CG = 0.1;
global CT = 0.1;
global GT = 0.3;

t = 1;
gtrmatrix = [*,AC*t,t,AT*t}{AC*t,*,CG*t,CT*t}{t,CG*t,*,GT*t}{AT*t,CT*t,GT*t,*];
HarvestFrequencies (vectorOfFrequencies,filteredData,1,1,1);
Model GTRmodel = (gtrmatrix,vectorOfFrequencies, 1);
Tree L_1 = ""+dummy;

/* Get scaling factor from GetString or manually*/
/*
GetString ( ModelString, GTRmodel, -1);
commandString = "";
commandString * 128;
commandString = "global scalingB := " + ModelString + ";";
ExecuteCommands ( commandString );
commandString * 0;
*/

/* manual */

global scalingB := 2*(vectorOfFrequencies[0):

*vectorOfFrequencies[1]*AC+vectorOfFrequencies[0]*vectorOfFrequencies[2]+           vectorOfFrequencies[0]*vectorOfFrequencies[3]*AT+vectorOfFrequencies[1]*vectorOfFrequencies[2]*CG+           vectorOfFrequencies[1]*vectorOfFrequencies[3]*CT+vectorOfFrequencies[2]*vectorOfFrequencies[3]*GT); /* Apply scaling factor */           for (k = 0; k < Columns (branchNames)-1; k=k+1)   {     ExecuteCommands ("L_1."+branchNames[k]+".t:="+branchLengthsNow[k]+"/scalingB;");   }   LikelihoodFunction lf = (filteredData,L_1);   Optimize (res,lf);   LIKELIHOOD_FUNCTION_OUTPUT=5;   fprintf (stdout, "\n",lf, "\n" );   branchLengthsAfter = BranchLength (L_1,-1); for ( k = 0; k < Columns ( branchNames) -1; k = k + 1 ) {      thisbranch = branchNames [ k ];      fprintf ( stdout, thisbranch, ": ", branchLengthsNow[ k ], " ", branchLengthsAfter [ k ], "\n" ); } fprintf (stdout, "\nBranch scaling factor computed to be ", scalingB,"\n\n");

Dear Wayne,

The GetString function actually returns the _expression_ for the branch length, including all parameters and such. The easiest way to strip out 't' is to do a regular expression search and replace in the returned string (using the '^' operator).

I'll work on a more portable way to do this...

Cheers,
Sergei

replaced t with 1 and now get convergence of manual scale factor with GetString. works well. thanks ./w

Dear Wayne,

Cool. What I will do for more portability is hook some analytical polynomial routines in HyPhy (a remnant of an old and never finished project) that will simplify the branch length expression automatically (if it is in fact a polynomial) and provide convenient factorizations for all parameters in which the expression is linear.

Cheers,
Sergei