Nick1979
YaBB Newbies
Offline
Curious HyPhy user
Posts: 6
Houston
Gender:

Hi Sergei,
I have two aligned sequences and i want to test whether the three codon sites evolve at equal rates.
So far I have the code below. Now I want to put a constraint in which sites 2 and 3 have the same branch length as site 1 and then compare the two optimized likelihood functions with an LRT. I will greatly appreciate if you help me with the coding ? Thank you
DataSet myData = ReadDataFile ("C:/Users/bio/Desktop/example");
/* Create 3 filters: one for each of the three codon positions */
DataSetFilter myFilter1 = CreateFilter (myData,1,"<100>","0,1,3"); DataSetFilter myFilter2 = CreateFilter (myData,1,"<010>","0,1,3"); DataSetFilter myFilter3 = CreateFilter (myData,1,"<001>","0,1,3");
HarvestFrequencies (obsFreqs1, myFilter1, 1, 1, 1); HarvestFrequencies (obsFreqs2, myFilter2, 1, 1, 1); HarvestFrequencies (obsFreqs3, myFilter3, 1, 1, 1);
/* Define an HKY85 model using each of the 4 frequency sets */ /* Let R be the ts/tv ratio, which we will eventually allow to be shared by all three of the positions */ global R;
HKY85RateMatrix = {{*,b,R*b,b}{b,*,b,R*b}{R*b,b,*,b}{b,R*b,b,*}}; Model HKY851 = (HKY85RateMatrix, obsFreqs1); Tree myTree1 = (a,b); Model HKY852 = (HKY85RateMatrix, obsFreqs2); Tree myTree2 = (a,b); Model HKY853 = (HKY85RateMatrix, obsFreqs3); Tree myTree3 = (a,b);
/* Analyze data using 3 partitions: each partition has separate rates and base frequencies, but the ts/tv ratio (R) is shared by all three. Also, each branch has the same ts/tv ratio. */
LikelihoodFunction theSplitLikFun = (myFilter1,myTree1,myFilter2,myTree2,myFilter3,myTree3);
Optimize (paramValues, theSplitLikFun); lnLik1 = paramValues[1][0]; npar1 = paramValues[1][1]+9;
