HyPhy message board - Modification of site specific rate estimation

	Welcome, Guest. Please Login

Home

Help

HyPhy message board › HYPHY Package › HyPhy feedback › Modification of site specific rate estimation

(Moderators: Sergei, Simon)

‹ Previous Topic | Next Topic ›

Pages: 1 2 3 4

Send Topic

Modification of site specific rate estimation (Read 20976 times)

Sergei

YaBB Administrator

Offline

Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male

Re: Modification of site specific rate estimation
Reply #30 - Mar 18^th, 2008 at 10:56pm

Dear Franscesc,

There have been so many revisions of this code:) Seems like there was some deprecated code that was messing up the units.

HyPhy will automatically unroot trees by default (and 'lose' a branch length); which is why you were seeing incorrect tree lengths reported. I added a line (ACCEPT_ROOTED_TREES) that forces HyPhy to use rooted trees directly.

I also posted a revised siteRates.bf script (Multimedia File Viewing and Clickable Links are available for Registered Members only!! You need to

Login) which should output rates in terms of units that you expect.
Please note that this script will only work with K2P (kappa=2), otherwise you have to compute another scaling factor that relates the internal model time parameter with the branch length (for K2P, kappa =2 this factor is 1).

Your our example now reports (note that I also added E[subs/site] to the output for your reference):

Code:

Reference tree length (units time): 20.62
Mean rate: 0.125607 substitutions/site/unit time
Mean substitutions/site: 2.59002
Site    1 Rate =  0.0548 E[Subs/site] =  1.1296 Log(L) -4.8524
Site    2 Rate =  0.0000 E[Subs/site] =  0.0000 Log(L) -1.3863
Site    3 Rate =  0.5083 E[Subs/site] = 10.4807 Log(L) -10.9215

Cheers,
Sergei

Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego

WWW

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #31 - Mar 19^th, 2008 at 3:02pm

Hi Sergei,
Thanks!
Some comments, let me know if I am wrong.
rateX is calculated from the tree using this line and it is the mean rate: substitutions/site/unit time.
Code:

ReplicateConstraint ("this1.?.?:=rateX*this2.?.?__",givenTree,givenTree);
LikelihoodFunction lf = (filteredData,givenTree);
Optimize (resNull,lf);

Is this double givenTree correct?

With this batch file, now, we multiply siteRate * rateX to get the number of substitution for that site and then we divided by the treeLength to get the rate.
i am confused about all this proces and the unist of siteRate which is what we optimize, right?

Some test: I repeated the same test I did before and everything looks ok. Except, one stupid thing. If I asked to show more decimal figures in the rates we can see slightly differences when we analyze 1 site or the whole alignment.
ex:
Quote:

Whole align. Site 1 Rate = 0.0547825 E[Subs/site] = 1.1296142 Log(L) -4.8524203
Single Site 1 Rate = 0.0216822 E[Subs/site] = 0.4470872 Log(L) -4.9924007
Whole align. Site 66 Rate = 3.7462178 E[Subs/site] = 77.2470103 Log(L) -15.2019927
Single Site 66 Rate = 3.7465842 E[Subs/site] = 77.2545666 Log(L) -15.2019927

There another question. I though if I change the model, especifying a new model and fixing all paramaters I could repeat the analisys.
Ex: I change the model by this one:
Code:

noequalFreqs = {{0.2897}{0.2316}{0.2425}{0.2362}};
GTRRateMatrix = {{*,b*3.1462,b*4.1238,b*1.1798}{b*3.1462,*,b*1.7147,b*4.8779}{b*4.1238,b*1.7147,*,b}{b*1.1798,b*1.7147,b,*}};
Model GTR = (GTRRateMatrix, noequalFreqs);
UseModel(GTR);

I am getting the following results:
Quote:

Reference tree length (units time): 20.62
Mean rate: 0.0706249 substitutions/site/unit time
Mean substitutions/site: 1.45629
Site 1 Rate = 0.0216822 E[Subs/site] = 0.4470872 Log(L) -4.9924007
Site 66 Rate = 75.3380671 E[Subs/site] = 1553.4709443 Log(L) -14.0156847

So, rateX changed, then it is not calculated from the tree. Is it calculated from the tree and the model?
I was expecting see a value higher than 1 for the site 1.

thanks,
Francesc

IP Logged

Sergei

YaBB Administrator

Offline

Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male

Re: Modification of site specific rate estimation
Reply #32 - Mar 19^th, 2008 at 4:00pm

Dear Francesc,

rateX is the 'mean' rate measures in substitutions/site/unit time.
siteRate is parameterized as a multiplicative factor of this mean rate. Could also be estimated directly, but it's just a matter of convenience.

Code:

ReplicateConstraint ("this1.?.?:=rateX*this2.?.?__",givenTree,givenTree);

is correct, because the __ in the constraint specs substitutes actual values from the chronogram (which are stored in givenTree to begin with).

You may be getting different results because OPTIMIZATION_PRECISION is too low. For a single site optimization you may need to set it to a small value in your script (e.g. OPTIMIZATION_PRECISION = 0.000001). It is set to 0.001 by default - that precision is on LogL values only (could be less on parameter values).

When you change the model to GTR, as I noted before, you need to compute another scaling factor, this one which maps 'b' to expected substitutions per site. You can do it like this (assuming 'b' is the branch length parameter and everything else in the model is constant)

Code:

noequalFreqs = {{0.2897}{0.2316}{0.2425}{0.2362}};
GTRRateMatrix = {{*,b*3.1462,b*4.1238,b*1.1798}{b*3.1462,*,b*1.7147,b*4.8779}{b*4.1238,b*1.7147,*,b}{b*1.1798,b*1.7147,b,*}};
Model GTR = (GTRRateMatrix, noequalFreqs);
GetString (bToSubs,GTR,-1); /* this returns an formula for the branch length in terms of b */
b=1;
ExecuteCommands ("bToSubs = " + bToSubs);
fprintf (stdout, "Scaling factor = ", bToSubs,"\n");

Then you need to multiply rateX by bToSubs. In this case bToSubs is about 1.83

HTH,
Sergei

Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego

WWW

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #33 - Mar 20^th, 2008 at 6:44am

Thanks for illuminating me!

I like to have the freedom to change the model, that was the reason I was testing GTR also.
So, if I want to try different models (hardcoded in the bach file), do I have to create this new scaling factor to multiply rateX? for all models except K2P?

So, the difference between the new siteRates.bf and the original one is that, in the new version, we have hardcoded the model and its parameters and, therefore, there is no estimation of the parameters.
In the siteRate.bf original, is it taking the tree and the topology as inputted or is is estimating branch length?

I just saw that when we use the original (like the new version) the rates printed in the log file are not equal to the ones printed in the window with the graphic.
Which ones are the correct in the original bf?

using the original template batch file siteRates.bf (JC69 model, fixed rates):
Log:
Quote:

Site 1 Rate = 0.4691 Log(L) -5.5499
Site 2 Rate = 0.0000 Log(L) -1.3863
Site 3 Rate = 3.9734 Log(L) -11.4171

window:
Quote:

rate Log
0.248794 -5.54995
0 -1.38629
2.10739 -11.4171

but, still, shouldn't we expect values for the 1st site of 0.05... (~1 substitution/20.62(tree length)).

thanks,
Francesc

« Last Edit: Mar 20^th, 2008 at 10:17am by Francesc »

IP Logged

Sergei

YaBB Administrator

Offline

Datamonkeys are forever...

Posts: 1658
UCSD
Gender: male

Re: Modification of site specific rate estimation
Reply #34 - Mar 20^th, 2008 at 10:05am

Dear Francesc,

You can use the same GetString trick to estimate the factor for other models (assuming they have 'b' as the branch length parameter).

The original siteRates.bf writes out rates relative to the mean (computed as the sum of site rates) to the log file (i.e. a value of 0.25 means that the rate is 25% of the mean rate for the alignment). This is to ensure that the mean rate in the output log is '1'.

What is printed to the console is the rate relative to the 'global' rate (estimated from the entire tree), i.e. you need to multiply the two to get absolute values.

Keep in mind that the original siteRates.bf was only intended to estimates site rates relative to each other Absolute values should not be interpreted directly, really, in the original version.

Cheers,
Sergei

Associate Professor
Division of Infectious Diseases
Division of Biomedical Informatics
School of Medicine
University of California San Diego

WWW

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #35 - Mar 31^st, 2008 at 9:43am

So, just to be sure I am doing the analyzes properly.... cos I am still confused when to use the rateX trick.
Don't I have the specify the trick always?

I don't have to use the trick if I specified the JC model

Code:

equalFreqs = {{0.25}{0.25}{0.25}{0.25}};
JCRateMatrix = {{*,b,b,b}{b,*,b,b}{b,b,*,b}{b,b,b,*}};
Model JC = (JCRateMatrix, equalFreqs);
UseModel(JC);

neither for K2P (kappa=2)

Code:

equalFreqs = {{0.25}{0.25}{0.25}{0.25}};
K2Pkappa2RateMatrix = {{*,b,b*2,b}{b,*,b,b*2}{b*2,b,*,b}{b,b*2,b,*}};
Model K2P = (K2Pkappa2RateMatrix, equalFreqs);
UseModel(K2P);

but I have to specified for others?
Code:

/*GTR*/
noequalFreqs = {{0.2897}{0.2316}{0.2425}{0.2362}};
GTRRateMatrix = {{*,b*3.1462,b*4.1238,b*1.1798}{b*3.1462,*,b*1.7147,b*4.8779}{b*4.1238,b*1.7147,*,b}{b*1.1798,b*1.7147,b,*}};
Model GTR = (GTRRateMatrix, noequalFreqs);
GetString (bToSubs,GTR,-1); /* this returns an formula for the branch length in terms of b */
b=1;
ExecuteCommands ("bToSubs = " + bToSubs);
fprintf (stdout, "Scaling factor = ", bToSubs,"\n");
.....more code......
rateX = rateX*bToSubs;


/*TrN*/
noequalFreqs = {{0.35374827}{0.16048399}{0.27201931}{0.21374843}};
TrNRateMatrix = {{*,b,b*2.21248665,b}{b,*,b,b*9.96807480}{b*2.21248665,b,*,b}{b,b*9.96807480,b,*}};
Model TrN = (TrNRateMatrix, noequalFreqs);
GetString (bToSubs,TrN,-1); /* this returns an formula for the branch length in terms of b */
b=1;
ExecuteCommands ("bToSubs = " + bToSubs);
fprintf (stdout, "Scaling factor = ", bToSubs,"\n");
.....more code......
rateX = rateX*bToSubs;

IP Logged

Art Poon

Global Moderator

Offline

Feed your monkey!

Posts: 0

Re: Modification of site specific rate estimation
Reply #36 - Mar 31^st, 2008 at 10:53am

Hi Francesc,

I'm covering for Sergei while he's in South Africa. I'm still trying to follow all that you guys have been discussing in this extensive thread, but as far as I understand there's no need to rescale branch lengths for the simpler models such as Kimura 2-parameter with kappa = 2 because all nucleotide states are evolving at the same rate (given uniform frequency distribution). There's nothing stopping you from computing rateX for JC or K2P(kappa=2), but you'll simply end up with rateX = 1. Wink

Best,
- Art.

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #37 - Apr 3^rd, 2008 at 8:24am

Thanks for your help.
So, given uniform frequency distribution of nucleotides calculation of bToSubs is not necessary.

I just realized about one problem with the file formats. I don't know if it is the right place to talk about it but it is related to this batchfile and my result
I am working with DNA sequences in nexus format and I am getting different results if I include the 'symbols' Format subcommand or not.
Is that a problem of hyphy reading the nexus format?
e.g.
With this file I am getting different results than
Code:

#NEXUS
begin data;
dimensions ntax=12 nchar=606;
format interleave datatype=dna   gap=-;

with this other file (the rest of the file is identical).
Code:

#NEXUS
begin data;
dimensions ntax=12 nchar=606;
format interleave datatype=dna   gap=- symbols="TCAG";

I am guessing that the correct values are without the symbol table....

Thanks a lot.

UPDATE:
If I use fasta format I am getting the same results as without symbol list.
If I use the default siteRates.bf, hyphy daesn't accept the nexus file with symbol list.
However, the error message doesn't specify what the problem is.

IP Logged

Art Poon

Global Moderator

Offline

Feed your monkey!

Posts: 0

Re: Modification of site specific rate estimation
Reply #38 - Apr 3^rd, 2008 at 10:43am

Dear Francesc,

The first thing that comes to mind is that the alphabet "ACGT" may be too restrictive. Are there any ambiguous base calls in your data?, e.g. N, R, Y, et cetera... Specifying symbols in a NEXUS file may influence how HyPhy is interpreting your sequences when this is the case. I would simply drop the symbols statement.

Best,
- Art

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #39 - Apr 3^rd, 2008 at 5:05pm

Hi Art,
yep, no problem with that.
I also wanted to note that there is no error message when opening a nexus file with the symbols list on it but the error message appears in other batchfiles.

I tested different models and it looks like it works fine.
Now, I would like to try to use aa models.
How should I hardcode the model? Using one of the pre-defined in hyphy (.mdl) would be fine.

I guess I should keep rateX as scaling factor (but not btoSubs of course).

Thanks again for your time and help.

IP Logged

Art Poon

Global Moderator

Offline

Feed your monkey!

Posts: 0

Re: Modification of site specific rate estimation
Reply #40 - Apr 3^rd, 2008 at 5:27pm

Hi Francesc,

Could you e-mail me your NEXUS file so that I can try to replicate the error on my machine? siteRates.bf uses ReadDataFile() just like any other batch file, so it's not obvious to me what could be causing the problem. I'd like to try breaking the batch file myself.

To use AA substitution models, I would convert your data file into AA sequences first and then call one of the template models in HyPhy.

Cheers,
- Art.

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #41 - Apr 4^th, 2008 at 9:56am

Hi,
Here you have the batchfile with a GTR model and the nexus file:
Multimedia File Viewing and Clickable Links are available for Registered Members only!! You need to

Login

For me it is not a problem, but maybe you want to know why these differences.
SiteRates.bf + nexus with symbols = message error
SiteRates.bf + nexus without symbols = works
SiteRates_GTR.bf + nexus with symbols = works
SiteRates_GTR.bf + nexus without symbols = works but different rates
I tried also K2P and I have the same results as with GTR but with JC I am getting the same rates with or without symbols.

Wink

yes, I have aa and nucleoides sequences. No problem with that.
My question was about how to call or hardcode a specific hyphy template aa model in the batchfile.
So, I want to change the nucleotide model (e.g. :
Code:

noequalFreqs = {{0.23632667}{0.29707150}{0.26400159}{0.20260024}};
GTRRateMatrix = {{*,b,b*4.05558988,b*1.35576178}{b,*,b*1.35576178,b*4.05558988}{b*4.05558988,b*1.35576178,*,b}{b*1.35576178,b*4.05558988,b,*}};
Model GTR = (GTRRateMatrix, noequalFreqs);
GetString (bToSubs,GTR,-1); /* this returns an formula for the branch length in terms of b */
b=1;
ExecuteCommands ("bToSubs = " + bToSubs);

)

for an aa model in the batch file.

Thanks!

IP Logged

Art Poon

Global Moderator

Offline

Feed your monkey!

Posts: 0

Re: Modification of site specific rate estimation
Reply #42 - Apr 4^th, 2008 at 12:05pm

Hi Francesc,

I wasn't able to download the batchfile from your link, so let's just go ahead with your amino acid problem.
As I understand it, you want to use a template AA model with hard-coded parameter values, as would be the case for an empirical rate matrix. Try the following code:

Code:

DataSet			ds = ReadDataFile ("/Applications/HyPhy/data/p51.aa");
DataSetFilter	dsf = CreateFilter (ds, 1);

SelectTemplateModel (dsf);
/* creates vector 'equalFreqs' which stores observed AA frequencies */

GetString (bToSubs,DayhoffModel,-1); /* this returns an formula for the branch length in terms of b */
t=1;
ExecuteCommands ("bToSubs = " + bToSubs);
fprintf (stdout, "bToSubs = ", bToSubs);

Note that this assumes that you know beforehand the name of the model that you're going to select, in this case "DayhoffModel", and that branch lengths use the local variable 't' rather than 'b', as in the previous example.

Cheers,
- Art.

IP Logged

Francesc

YaBB Newbies

Offline

I love YaBB 1G - SP1!

Posts: 32

Re: Modification of site specific rate estimation
Reply #43 - Apr 4^th, 2008 at 1:53pm

Here they are again:
Multimedia File Viewing and Clickable Links are available for Registered Members only!! You need to

Login

Concerning to the aa models...
I am doing this exercises, not only cos I have to analyze some data, but also cos I want to pipeline hyphy with other programs I am preparing.

I tried the following using as a base (siteRates_GTR):
Code:

SetDialogPrompt ("Please specify a amino-acid data file:");

DataSet ds = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter filteredData = CreateFilter (ds,1);

fprintf (stdout,"\n______________READ THE FOLLOWING DATA______________\n",ds);

categoriesUsed = 0;

SelectTemplateModel (filteredData);
/* creates vector 'equalFreqs' which stores observed AA frequencies */

GetString (bToSubs,DayhoffModel,-1); /* this returns an formula for the branch length in terms of b */
t=1;
ExecuteCommands ("bToSubs = " + bToSubs);
fprintf (stdout, "bToSubs = ", bToSubs);

It asks for aa file but it also asks for an aa model. ??
't' doesn't appear in the code again. So, is it doing something?

Assuming I want to use a model without using empirical base frequencies.
Code:

SetDialogPrompt ("Please specify amino-acid data file:");

DataSet ds = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter filteredData = CreateFilter (ds,1);

fprintf (stdout,"\n______________READ THE FOLLOWING DATA______________\n",ds);

categoriesUsed = 0;

UseModel(DayhoffModel);

This doesn't work right?
I removed bToSubs cos it wouldn't be necessary anymore. We only used when aa frequencies are not uniform, right?

Sorry, I guess it is not easy to arrived in middle of the topic.

Thanks,
Francesc

IP Logged

Art Poon

Global Moderator

Offline

Feed your monkey!

Posts: 0

Re: Modification of site specific rate estimation
Reply #44 - Apr 4^th, 2008 at 2:11pm

Hi Francesc,

Try this:

Code:

DataSet			ds = ReadDataFile ("/Applications/HyPhy/data/p51.aa");
DataSetFilter	dsf = CreateFilter (ds, 1);

equalFreqs = {20,1}+0.05;

myMatrix = {20,20};

/* fill these with your own values, or, if many cells are constrained to the same value
    then you can be clever and write a loop statement over those cells  */

_numericRateMatrix =
{
{0,	0.167653,	4.43521,	5.56325,	0.597923,	1.8685,	0.005,	0.005,	0.592784,	0.16024,	0.005,	0.617509,	1.00981,	0.005,	0.0744808,	8.5942,	24.1422,	24.8094,	0.005,	0.005}
{0.167653,	0,	0.005,	0.005,	0.362959,	0.0489798,	0.005,	0.005,	0.005,	0.005,	0.005,	0.0604932,	0.005,	0.005,	2.86364,	1.12195,	0.005,	0.005,	5.49894,	8.34835}
{4.43521,	0.005,	0,	12.1233,	0.005,	10.3969,	2.31779,	0.145124,	0.894313,	0.005,	0.005,	29.4087,	0.005,	0.005,	0.0674539,	0.427881,	0.630395,	2.91786,	0.005,	2.28154}
{5.56325,	0.005,	12.1233,	0,	0.005,	14.7801,	0.005,	0.0390512,	23.9626,	0.129839,	0.005,	0.201526,	0.005,	3.20656,	0.0251632,	0.005,	0.458743,	2.19952,	0.005,	0.005}
{0.597923,	0.362959,	0.005,	0.005,	0,	0.005,	0.005,	1.48288,	0.005,	7.48781,	0.005,	0.005,	0.0342252,	0.005,	0.005,	4.27939,	0.114512,	2.28,	0.005,	4.12728}
{1.8685,	0.0489798,	10.3969,	14.7801,	0.005,	0,	0.005,	0.005,	0.279425,	0.0489798,	0.0489798,	0.0604932,	0.005,	0.0604932,	13.4379,	6.27966,	0.0489798,	2.79622,	2.8258,	0.005}
{0.005,	0.005,	2.31779,	0.005,	0.005,	0.005,	0,	0.005,	0.22406,	1.76382,	0.005,	8.59876,	13.9444,	18.5465,	6.84405,	0.725157,	0.95956,	0.827479,	0.005,	47.4889}
{0.005,	0.005,	0.145124,	0.0390512,	1.48288,	0.005,	0.005,	0,	0.817481,	9.10246,	17.3064,	0.987028,	0.005,	0.0342252,	1.34069,	0.740091,	9.36345,	24.8231,	0.005,	0.114512}
{0.592784,	0.005,	0.894313,	23.9626,	0.005,	0.279425,	0.22406,	0.817481,	0,	0.005,	4.09564,	10.6655,	0.111928,	13.0705,	39.8897,	0.005,	4.04802,	0.128065,	0.005,	0.005}
{0.16024,	0.005,	0.005,	0.129839,	7.48781,	0.0489798,	1.76382,	9.10246,	0.005,	0,	11.3839,	0.005,	9.83095,	2.89048,	0.586757,	6.14396,	0.005,	2.95344,	1.37031,	0.005}
{0.005,	0.005,	0.005,	0.005,	0.005,	0.0489798,	0.005,	17.3064,	4.09564,	11.3839,	0,	0.201526,	0.005,	0.005,	3.28652,	0.392575,	7.41313,	14.7683,	0.005,	0.579198}
{0.617509,	0.0604932,	29.4087,	0.201526,	0.005,	0.0604932,	8.59876,	0.987028,	10.6655,	0.005,	0.201526,	0,	0.344848,	0.342068,	0.16024,	14.5699,	4.54206,	0.0744808,	0.005,	5.06475}
{1.00981,	0.005,	0.005,	0.005,	0.0342252,	0.005,	13.9444,	0.005,	0.111928,	9.83095,	0.005,	0.344848,	0,	3.04502,	0.404723,	14.249,	4.33701,	0.005,	0.005,	0.005}
{0.005,	0.005,	0.005,	3.20656,	0.005,	0.0604932,	18.5465,	0.0342252,	13.0705,	2.89048,	0.005,	0.342068,	3.04502,	0,	10.6746,	0.16024,	0.203091,	0.005,	0.0443298,	0.005}
{0.0744808,	2.86364,	0.0674539,	0.0251632,	0.005,	13.4379,	6.84405,	1.34069,	39.8897,	0.586757,	3.28652,	0.16024,	0.404723,	10.6746,	0,	8.35024,	0.928203,	0.279425,	5.96564,	0.005}
{8.5942,	1.12195,	0.427881,	0.005,	4.27939,	6.27966,	0.725157,	0.740091,	0.005,	6.14396,	0.392575,	14.5699,	14.249,	0.16024,	8.35024,	0,	6.34079,	0.862637,	1.10156,	0.933142}
{24.1422,	0.005,	0.630395,	0.458743,	0.114512,	0.0489798,	0.95956,	9.36345,	4.04802,	0.005,	7.41313,	4.54206,	4.33701,	0.203091,	0.928203,	6.34079,	0,	0.005,	0.005,	0.490608}
{24.8094,	0.005,	2.91786,	2.19952,	2.28,	2.79622,	0.827479,	24.8231,	0.128065,	2.95344,	14.7683,	0.0744808,	0.005,	0.005,	0.279425,	0.862637,	0.005,	0,	0.005,	1.35482}
{0.005,	5.49894,	0.005,	0.005,	0.005,	2.8258,	0.005,	0.005,	0.005,	1.37031,	0.005,	0.005,	0.005,	0.0443298,	5.96564,	1.10156,	0.005,	0.005,	0,	0.005}
{0.005,	8.34835,	2.28154,	0.005,	4.12728,	0.005,	47.4889,	0.114512,	0.005,	0.005,	0.579198,	5.06475,	0.005,	0.005,	0.005,	0.933142,	0.490608,	1.35482,	0.005,	0}
};


/*--------------------------------------------------------------------------------------------*/

function PopulateModelMatrix (ModelMatrixName&, EFV)
{
	ModelMatrixName = {20,20};
	if (categoriesUsed)
	{
		for (ri = 0; ri < 20; ri = ri+1)
		{
			for (ci = ri+1; ci < 20; ci = ci+1)
			{
				ModelMatrixName[ri][ci] := _numericRateMatrix__[ri__][ci__] * c * t;
				ModelMatrixName[ci][ri] := _numericRateMatrix__[ri__][ci__] * c * t;
			}
		}
	}
	else
	{
		for (ri = 0; ri < 20; ri = ri+1)
		{
			for (ci = ri+1; ci < 20; ci = ci+1)
			{
				ModelMatrixName[ri][ci] := _numericRateMatrix__[ri__][ci__] * t;
				ModelMatrixName[ci][ri] := _numericRateMatrix__[ri__][ci__] * t;
			}
		}
	}
	return 1;
}

/*--------------------------------------------------------------------------------------------*/


MULTIPLY_BY_FREQS = PopulateModelMatrix ("myMatrix",equalFreqs);

Model myModel = (myMatrix, equalFreqs, MULTIPLY_BY_FREQS);

FREQUENCY_SENSITIVE = 0;
UseModel (myModel);

Cheers,
- Art.

IP Logged

Pages: 1 2 3 4

Send Topic

‹ Previous Topic | Next Topic ›

« Home

‹ Board

Top of this page