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):
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. |