Dear all,

It took me some time figuring out how to build a customized batch file for HyPhy with a REV matrix, plus a proportion of invariable sites and the other ones following gamma... I hope this can be useful to someone else as well (tested with the same tree and same parameters estimated from PAUP*, same results were obtained):

-----------------------------------------------------------------
/* Basefreqs, 5 transition types, proportion of invariable sites
  and rate heterogeneity are global. */

#include "displayfunction.bf";

myTopology = "some_tree"

DataSet myData  = ReadDataFile("some_data");
DataSetFilter myFilter = CreateFilter(myData,1);
HarvestFrequencies(obsFreqs, myFilter, 1, 1, 1);

/* "invariant sites" rate and "gamma shape (alpha) for rate heterogeneity" categories */
global Pinvar;
global shape;

/* Five individual global transition rates of each type (the other one is free to vary on each branch). */
global r1;
global r2;
global r3;
global r4;
global r5;

/* Discrete probability distribution used to estimate # of invariant sites:
  2 classes, 1 with invariant sites (0), the other with variable ones (1).
  The value of "Pinvar" shall be within the interval [0,1]
  Estimate 'Pinvar' */
category invClass = (2,[Pinvar}{1-Pinvar],MEAN, ,[0}{1], 0,1);

/* Gamma distributed site rate heterogenity:
  4 classes, gamma prob. distribution of rates.
  Estimate 'shape' (= alpha) */
category rateCat =(4, EQUAL, MEAN, GammaDist(_x_,shape,shape), CGammaDist(_x_,shape,shape),0,1e25);

/* REV+I+G matrix. (I let the "C-T" rate to freely vary, but could be any other as well) */
REV_I_G_rateMatrix =
[         *          ,r1*m*invClass*rateCat,r2*m*invClass*rateCat,r3*m*invClass*rateCat}
{r1*m*invClass*rateCat,         *           ,r4*m*invClass*rateCat,   m*invClass*rateCat}
{r2*m*invClass*rateCat,r4*m*invClass*rateCat,         *           ,r5*m*invClass*rateCat}
{r3*m*invClass*rateCat,   m*invClass*rateCat,r5*m*invClass*rateCat,         *           ];

Model REV_I_G = (REV_I_G_rateMatrix, obsFreqs);

UseModel(REV_I_G);
Tree myTree = myTopology;
LikelihoodFunction theLikFun = (myFilter, myTree);
Optimize(results,theLikFun);
fprintf (stdout,theLikFun);

/* Log 'Pinvar' to screen. */
GetInformation (catInfo,invClass);
catInfo = echoinvClass(catInfo);

/* Log 'shape' to screen. */
GetInformation (catInfo,rateCat);
catInfo = echoinvClass(catInfo);

---------------------------------------------------------------

Dear Jose,

Thanks! I haven't thought of this way to parameterize the G+I (using multiplicative category variables), it's actually a neat trick.

For the information of other users: 'displayfunction.bf' is located in Examples/Categories.

Cheers,
Sergei

hi,

I tried to use this trick in my own code and I got the following error message:

ERROR: only one category variable can be handled by this code

did I miss something or did somethind change in the way hyphy handle multiple category variables?

thanks,

Olivier

Dear Olivier,

I think you may need to redefine the 'echoinvClass' function (in displayfunction.bf) to do this. It's not difficult to do, because all you have to do is add a '0' rate with weight Pinvar to the display, and then multiply all gamma weights by (1-Pinvar).

Cheers,
Sergei

Thanks for your email Sergei. Sorry, I looked more carefully and I located the problem. Actually I tried to plug in this G+I code into a modified Plato.bf. But the plato.bf does not handle multiple category variable (?). Is there an other way to programatically define REV+G+I?

Thanks,

Olivier

Dear Olivier,

You can define the G+I category variable (resp bins in Gamma) thus:


Code (]
     global alpha   = .5;
     alpha:>0.01;alpha:<100;

     global beta = 1;
     beta:>0.01;
     beta:<200;

     global P;
     P :< 1-(1e-10);
     P = 1/(resp+1);
     
     beta := (1-P)*alpha;
     
     
     catFreqMatrix = {resp+1,1};
     catFreqMatrix [0):

:= P;            for (h=1; h<=resp;h=h+1)      {            catFreqMatrix [h] := (1-P)/resp__;      }                                      category c = (resp+1, catFreqMatrix, MEAN,                              (1-P)*GammaDist(_x_,alpha,beta)*(_x_>0),                              (1-P)*CGammaDist(_x_,alpha,beta)*(_x_>0)+P,                              0,                              1e25,                              (1-P)*CGammaDist(_x_,alpha+1,beta)*(alpha/beta)*(_x_>0)                           );

Cheers,
Sergei

Hi Sergei, this could deserve a new thread but is nonetheless related to this one. How is the conditional site likelihood matrix constructed when using multiplicative category variables and subsequently calling ConstructCategoryMatrix. I have two category variables (3 and 2 intervals respectively) which gives a (3*2)*sites matrix but some values are repeating. For instance
{ { 4.904 e-05, ...}
  { 4.904 e-05, ...}
  {3.480 e-05, ....}
  {3.480 e-05, ...}
   etc
}
thanks ./w

Dear Wayne,

For multiple cat. variables, each column of 'ConstructCategoryMatrix' will be a flattened representation of a multi-dimensional matrix of all possible cat. var. values.

Each likelihood function has a ordering of category variables (usually in the order they were encountered during likelihood function construction), which can be retrieved at run time using

[code]
GetInformation (storageID, lfID);
[/code]

Suppose for example, that there are two category variables: c1 (M classes) and c2 (N classes) and c1 is the first in the orderings.

Then the i-th entry (0-based indexing) of each column vector from 'ConstructCategoryMatrix' will correspond to i div N class in c1 and i mod N class in c2.

In your case, repeating values suggest that the second category variable does not affect that particular site (for whatever reason).

Cheers,
Sergei

thanks Sergei, i suspected a problem with my 2nd cat variable. now confirmed. thanks ./w

Hi Sergei,
On a related note, how do I specify an HKY85, beta-gamma (5 categories) model in a batch file (or for that matter, any standard model + beta-gamma)? I am trying to implement a sliding window NT model fit for command-line submission on a cluster. If there is an example or documentation that I missed, please feel free to point me there. Thanks in advance.

Regards,
Kristen Dang

Dear Kristen,


kdang wrote on Oct 16^th, 2007 at 6:24pm:

On a related note, how do I specify an HKY85, beta-gamma (5 categories) model in a batch file (or for that matter, any standard model + beta-gamma)? I am trying to implement a sliding window NT model fit for command-line submission on a cluster. If there is an example or documentation that I missed, please feel free to point me there. Thanks in advance.

Do you want to embed the model definition in your file?
If not, just select

'Rate variation'
'Beta-Gamma'
5

when a standard analysis prompts you for model options.
You can define HKY85 with beta-gamma with this self contained code (replace < insert_your_filter_name_here> as appropriate):

[code]
       resp = 5; /* how many categories */
     global betaP = 1;
     global betaQ = 1;
     betaP:>0.05;betaP:<85;
     betaQ:>0.05;betaQ:<85;
     category pc = (resp-1, EQUAL, MEAN,
                             _x_^(betaP-1)*(1-_x_)^(betaQ-1)/Beta(betaP,betaQ), /* density */
                             IBeta(_x_,betaP,betaQ), /*CDF*/
                             0,                           /*left bound*/
                             1,                     /*right bound*/
                             IBeta(_x_,betaP+1,betaQ)*betaP/(betaP+betaQ)
                          );
     
     global alpha = .5;
     alpha:>0.01;alpha:<100;
     category c = (resp, pc, MEAN,
                             GammaDist(_x_,alpha,alpha),
                             CGammaDist(_x_,alpha,alpha),
                             0 ,
                             1e25,
                             CGammaDist(_x_,alpha+1,alpha)
                          );


     global kappa_inv = 1;
     HKY_Rate_Matrix = [*,kappa_inv*t*c,t*c,kappa_inv*t*c}{kappa_inv*t*c,*,kappa_inv*t*c,t*c}{t*c,kappa_inv*t*c,*,kappa_inv*t*c}{kappa_inv*t*c,t*c,kappa_inv*t*c,*];
     
   HarvestFrequencies (vectorOfFrequencies, <insert_your_filter_name_here> ,1,1,0);
   Model HKY85Model = (HKY_Rate_Matrix, vectorOfFrequencies, 1);

[/code]

Cheers,
Sergei

Yes, I do want to embed it in the file (I'm trying to remove the interactive feature so I can submit to the parallel queue on our cluster).

Thanks! :)

Dear Kristen,

You can also use 'ExecuteAFile' with appropriately constructed options.
This saves a lot of copy and paste and makes the code more modular.

See Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login

Cheers,
Sergei