Jose_Patane
YaBB Newbies
Offline
Posts: 26
Brazil
Gender:

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}{1Pinvar}},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 "CT" 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);

