Welcome, Guest. Please Login
 YaBB - Yet another Bulletin Board
 
 
  HomeHelpSearchLogin  
 
  HyPhy message board › HYPHY Package › HyPhy feedback › generic branch scaling function
(Moderators: Sergei, Simon)
‹ Previous Topic | Next Topic ›
Page Index Toggle Pages: 1
generic branch scaling function (Read 4133 times)
wayne
Global Moderator
*****
Offline


I love YaBB 1G - SP1!

Posts: 57
San Diego, CA
Gender: male
generic branch scaling function
Feb 26th, 2007 at 2:09am  
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
Back to top
 
Assistant Project Scientist&&Antiviral Research Center&&Department of Pathology&&University of California, San Diego
WWW WWW  
IP Logged
 
wayne
Global Moderator
*****
Offline


I love YaBB 1G - SP1!

Posts: 57
San Diego, CA
Gender: male
Re: generic branch scaling function
Reply #1 - Feb 26th, 2007 at 9:57am  
Oops, the "stringB * 0;" in the third line of code shouldn't be there. ./w
Back to top
 
Assistant Project Scientist&&Antiviral Research Center&&Department of Pathology&&University of California, San Diego
WWW WWW  
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: generic branch scaling function
Reply #2 - Feb 26th, 2007 at 5:32pm  
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
Back to top
 
Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: generic branch scaling function
Reply #3 - Mar 12th, 2007 at 5:02pm  
Dear Wayne,

As of today's build calling

Code:
GetString(stringID,model_id,-1)


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

Cheers,
Sergei
Back to top
 
Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
wayne
Global Moderator
*****
Offline


I love YaBB 1G - SP1!

Posts: 57
San Diego, CA
Gender: male
Re: generic branch scaling function
Reply #4 - Mar 15th, 2007 at 1:31am  
thanks ./w
Back to top
 
Assistant Project Scientist&&Antiviral Research Center&&Department of Pathology&&University of California, San Diego
WWW WWW  
IP Logged
 
wayne
Global Moderator
*****
Offline


I love YaBB 1G - SP1!

Posts: 57
San Diego, CA
Gender: male
Re: generic branch scaling function
Reply #5 - Mar 19th, 2007 at 1:53am  
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");
Back to top
 
Assistant Project Scientist&&Antiviral Research Center&&Department of Pathology&&University of California, San Diego
WWW WWW  
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: generic branch scaling function
Reply #6 - Mar 19th, 2007 at 1:58pm  
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
Back to top
 
Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
wayne
Global Moderator
*****
Offline


I love YaBB 1G - SP1!

Posts: 57
San Diego, CA
Gender: male
Re: generic branch scaling function
Reply #7 - Mar 19th, 2007 at 3:10pm  
replaced t with 1 and now get convergence of manual scale factor with GetString. works well. thanks ./w
Back to top
 
Assistant Project Scientist&&Antiviral Research Center&&Department of Pathology&&University of California, San Diego
WWW WWW  
IP Logged
 
Sergei
YaBB Administrator
*****
Offline


Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male
Re: generic branch scaling function
Reply #8 - Mar 19th, 2007 at 4:33pm  
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
Back to top
 
Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego
WWW WWW  
IP Logged
 
Page Index Toggle Pages: 1
‹ Previous Topic | Next Topic ›
 
  « Home ‹ Board Top of this page

HyPhy message board » Powered by YaBB 2.5.2!
YaBB Forum Software © 2000-2024. All Rights Reserved.

Valid RSS Valid XHTML Valid CSS Powered by Perl Source Forge