HyPhy message board
http://www.hyphy.org/cgi-bin/hyphy_forums/YaBB.pl
HYPHY Package >> HyPhy bugs >> fscanf problem
http://www.hyphy.org/cgi-bin/hyphy_forums/YaBB.pl?num=1171390223

Message started by scottdoniger on Feb 13th, 2007 at 10:10am

Title: fscanf problem
Post by scottdoniger on Feb 13th, 2007 at 10:10am
Hi there - I'm trying to read in nucleotide frequencies from a file (I want to specifically set these for the file, rather than using HarvestFreqs). So I have a file with 4 frequencies, e.g.,

.001
.498
.001
.498

then in the batch file I have:
fscanf(stdin, "Number", F1);
fscanf(stdin, "Number", F2);
fscanf(stdin, "Number", F3);
fscanf(stdin, "Number", F4);
partition3_Freqs = {4,1};
partition3_Freqs[0][0]:= F1;
partition3_Freqs[1][0]:= F2;
partition3_Freqs[2][0]:= F3;
partition3_Freqs[3][0]:= F4;

fprintf(stdout, "A = ",partition3_Freqs[0][0],"\n");
fprintf(stdout, "C = ",partition3_Freqs[1][0],"\n");
fprintf(stdout, "G = ",partition3_Freqs[2][0],"\n");
fprintf(stdout, "T = ",partition3_Freqs[3][0],"\n");

HYPHY correctly prints out the frequencies, so it's reading it from the file, but
for some reason the likelihood is actually twice the likelihood when I set the frequencies by hand. That is if I type partition_Freqs[0][0] := .001, I get a different result than if I read them from the file.

It seems like the difference is a result of the fact that when I read in the nuc. frequencies, HYPHY does not change the TVTS ratio from it's starting point. When I set the frequencies by hand, it estimates TVTS as 0.

Here's the rest of my batch file.
fscanf(stdin, "String", file);
DataSet b = ReadDataFile(file);
DataSetFilter partition3 = CreateFilter(b,1);
global partition3_Shared_TVTS = .209335;
partition3_HKY85={4,4};
partition3_HKY85[0][1]:=t*partition3_Shared_TVTS;
partition3_HKY85[0][2]:=t;
partition3_HKY85[0][3]:=t*partition3_Shared_TVTS;
partition3_HKY85[1][0]:=t*partition3_Shared_TVTS;
partition3_HKY85[1][2]:=t*partition3_Shared_TVTS;
partition3_HKY85[1][3]:=t;
partition3_HKY85[2][0]:=t;
partition3_HKY85[2][1]:=t*partition3_Shared_TVTS;
partition3_HKY85[2][3]:=t*partition3_Shared_TVTS;
partition3_HKY85[3][0]:=t*partition3_Shared_TVTS;
partition3_HKY85[3][1]:=t;
partition3_HKY85[3][2]:=t*partition3_Shared_TVTS;


fscanf(stdin, "Number", F1);
fscanf(stdin, "Number", F2);
fscanf(stdin, "Number", F3);
fscanf(stdin, "Number", F4);
partition3_Freqs = {4,1};
/*
partition3_Freqs[0][0]:= F1;
partition3_Freqs[1][0]:= F2;
partition3_Freqs[2][0]:= F3;
partition3_Freqs[3][0]:= F4;
*/

partition3_Freqs[0][0]:= .001;
partition3_Freqs[1][0]:= .498;
partition3_Freqs[2][0]:= .001;
partition3_Freqs[3][0]:= .498;

fprintf(stdout, "A = ",partition3_Freqs[0][0],"\n");
fprintf(stdout, "C = ",partition3_Freqs[1][0],"\n");
fprintf(stdout, "G = ",partition3_Freqs[2][0],"\n");
fprintf(stdout, "T = ",partition3_Freqs[3][0],"\n");

Model partition3_HKY85_model2=(partition3_HKY85, partition3_Freqs);
global r2 = 1;
UseModel (partition3_HKY85_model2);
Tree Noncoding2 =(((Scer,Spar)Node2,Smik)Node1,Skud,Sbay);
Noncoding2.Scer.t=0;
Noncoding2.Spar.t=0;
Noncoding2.Node2.t=0;
Noncoding2.Smik.t=0;
Noncoding2.Node1.t=0;
Noncoding2.Skud.t=0;
Noncoding2.Sbay.t=0;

LikelihoodFunction test_LF2 = (partition3, Noncoding2);

Optimize(res_test_LF2,test_LF2);
fprintf(stdout,"\n",test_LF2);

My sequence files looks like:
>Scer
CTTT
>Spa r
CTCT
>Smik
TTTT
>Skud
TTTT
>Sbay
TTTT

Title: Re: fscanf problem
Post by Sergei on Feb 13th, 2007 at 10:36am
Dear Scott,

When you use the assign operator (:=) to set the frequencies based on F1-F4 parameters, what is actually happening is that HyPhy will use F1,F2,F3,F4 as model parameters to be estimated, hence you are actually fitting a different model. If you print partition3_Freqs  after a call to Optimize, you will find that the values have changed. Try using the set value (=) operator instead:


Code (]
partition3_Freqs[0):

[0]= F1;
partition3_Freqs[1][0]= F2;
partition3_Freqs[2][0]= F3;
partition3_Freqs[3][0]= F4;


Cheers,
Sergei

Title: Re: fscanf problem
Post by scottdoniger on Feb 13th, 2007 at 10:55am
ah. ok. thanks.

Now I've got a question about = vs. :=.

If I want to force my likelihood function to use a particular value (e.g. a TVTS ratio or substitution rate) I say something like:

global partition3_Shared_TVTS := .209335;

or
Synonymous.Scer.t:=.649815;

Is that correct? Because I'm assigning these parameters a number rather than a variable, they are now fixed.

If I say Synonymous.Scer.t =.649815; The model will use .649815 as a starting point and then estimate t for this branch, right?

thanks again. HYPHY is great!

Scott

Title: Re: fscanf problem
Post by Sergei on Feb 13th, 2007 at 11:40am
Dear Scott,

Yes, you are correct. ':= number' will assign a (fixed value) constraint to a variable, whereas '= number' will simply set the initial value for optimization.

Cheers,
Sergei

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