jamescotton
YaBB Newbies
Offline
Posts: 28
London
Gender:
|
Hi Sergei - thanks for the previous suggestions.. I've made a little bit of progress with understanding HyPhy, but not enough, i'm afraid.. I've got this far, and have a problem in that the branch lengths stored in the tree object don't seem to be updated.. I get output like this:
[quote] branch scaling factor computed to be 0.161777
BLmod=0.044344 Branch Pe_occidentalis. Length before scaling: 0.00717 Length after scaling : 0.00717 BLmod=0.077582 Branch Pe_conspicillatus. Length before scaling: 0.01255 Length after scaling : 0.01255 BLmod=0.100006 Branch M_serrator. Length before scaling: 0.01618 Length after scaling : 0.01618 BLmod=0.0441 Branch S_sula. Length before scaling: 0.00713 Length after scaling : 0.00713 [/quote]
and so on for every branch - the BLmod values are correct, but don't seem to store in the tree properly. The same happens if i use '=' rather than ':=' for the assignment in the first "executecommands" statement in the loop. I don't understand why this is happening!.. The second problem is that all the branch lengths are estimated to be the same length in the output.. [quote] global AC=0.2932585265898854; global AT=0.1951842034375856; global CG=0.07364278346962116; global CT=1.274516987700602; global GT=0.03398007119957082; k=0.02; scalingB=0.05617924196392339; L_1.Pe_occidentalis.t:=branchLengthsNow[k]/scalingB; L_1.Pe_conspicillatus.t:=branchLengthsNow[k]/scalingB; L_1.M_serrator.t:=branchLengthsNow[k]/scalingB; L_1.S_sula.t:=branchLengthsNow[k]/scalingB; L_1.S_dactylatra.t:=branchLengthsNow[k]/scalingB; L_1.Node7.t:=branchLengthsNow[k]/scalingB; L_1.Node5.t:=branchLengthsNow[k]/scalingB; L_1.A_anhinga.t:=branchLengthsNow[k]/scalingB; L_1.A_novaehollandiae.t:=branchLengthsNow[k]/scalingB; L_1.A_rufa.t:=branchLengthsNow[k]/scalingB; L_1.Node12.t:=branchLengthsNow[k]/scalingB; L_1.Node10.t:=branchLengthsNow[k]/scalingB; L_1.Node4.t:=branchLengthsNow[k]/scalingB; L_1.P_carbo.t:=branchLengthsNow[k]/scalingB; L_1.P_varius.t:=branchLengthsNow[k]/scalingB; L_1.Node16.t:=branchLengthsNow[k]/scalingB; L_1.P_melanoleucos.t:=branchLengthsNow[k]/scalingB; L_1.Node15.t:=branchLengthsNow[k]/scalingB; L_1.Node3.t:=branchLengthsNow[k]/scalingB; Tree L_1=(Pe_occidentalis:0.0469995,Pe_conspicillatus:0.0469995,(((M_serrator:0.04699 95,(S_sula:0.0469995,S_dactylatra:0.0469995)Node7:0.0469995)Node5:0.0469995,(A_a nhinga:0.0469995,(A_novaehollandiae:0.0469995,A_rufa:0.0469995)Node12:0.0469995) Node10:0.0469995)Node4:0.0469995,((P_carbo:0.0469995,P_varius:0.0469995)Node16:0 .0469995,P_melanoleucos:0.0469995)Node15:0.0469995)Node3:0.0469995); [/quote]
I *think* I understand this - i'm constraining each branch length to the value of branchLengthsNow[k]/scalingB, and the program tries to optimise this parameter (scalingB) and so optimises all the branch lengths to be the same.. do I need to somehow "de-reference" these values to say "constrain the branch length to the current value of branchLengthNow[h]/scalingB rather than to be the value of these variables that can then be optimised? i've no idea how to do this, or even if this makes sense in HyPhy..
My program is pasted below - sorry for such a long post, but i'm finding it a rather steep learning curve just now!
Thanks, James [quote] DataSet ds = ReadDataFile ("Darterfile.nex"); DataSetFilter filteredData = CreateFilter (ds,1);
ACCEPT_BRANCH_LENGTHS = 1; ACCEPT_ROOTED_TREES = 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); /* NOT SURE IF THIS SHOULD BE 1 or 0 FOR THE LAST PARAMETER */ Tree L_1 = "(Pe_occidentalis:0.044344,Pe_conspicillatus:0.077582,(((M_serrator:0.100006,(S_ sula:0.044100,S_dactylatra:0.076275):0.103486):0.037020,(A_anhinga:0.197727,(A_n ovaehollandiae:0.106082,A_rufa:0.068100):0.071018):0.079408):0.026938,((P_carbo: 0.033352,P_varius:0.037379):0.034797,P_melanoleucos:0.101265):0.131502):0.223011 );";
scalingB = computeScalingFactorB (gtrmatrix, vectorOfFrequencies); fprintf (stdout, "\nBranch scaling factor computed to be ", scalingB,"\n\n");
branchNames = BranchName (L_1,-1); branchLengthsNow = BranchLength (L_1,-1);
for (k = 0; k < Columns (branchNames)-1; k=k+1) { ExecuteCommands ("L_1."+branchNames[k]+".t:=branchLengthsNow[k]/scalingB;"); fprintf (stdout,"BLmod=",branchLengthsNow[k]/scalingB,"\n"); fprintf (stdout, "Branch ", branchNames[k],".\n\t Length before scaling: ", Format (branchLengthsNow[k],10,5), "\n\t Length after scaling : ", Format (BranchLength (L_1,k),10,5), "\n"); }
LikelihoodFunction lf = (filteredData,L_1);
Optimize (res,lf); LIKELIHOOD_FUNCTION_OUTPUT=5; fprintf (stdout, "\n",lf);
function computeScalingFactorB (rateMatrix, baseFreqs) { B = 0; for (n1 = 0; n1 < 4; n1 = n1+1) { for (n2 = 0; n2 < 4; n2 = n2+1) { if (n2!=n1) { B = B + baseFreqs[n1]*baseFreqs[n2]*rateMatrix[n1][n2]; } } } return B; } [/quote]
|