1/* This file defines the transition matrix for the Muse-Gaut 94 model.
2   The file should be used as follows:
3
4   1) Read Data File and create datafilter filteredData
5   2) #include this file (or use SelectTemplateModel(filteredData);)
6   3) Define the tree
7   4) Proceed with the likelihood function using 'vectorOfFrequencies' as the vector to pass to the constructor.
8
9   This model has the following signature:
10   	#Short:MG94#
11   	#Desc:Muse-Gaut 94 codon model. Local or global parameters. Possible Gamma Variation.#
12   	#Dimension:*#
13    #DataType:codon#
14   	#FileName:MG94.mdl#
15
16   08/18/1999  by Sergei L. Kosakovsky Pond
17*/
18
19ModelMatrixDimension = 0;
20
21function GetBranchDNDS (shortName)
22{
23	sR  = "givenTree."+shortName+".synRate";
24	nsR = "givenTree."+shortName+".nonSynRate";
25	sR = valueGrab (sR);
26	nsR = valueGrab (nsR);
27	if (sR > 0.0)
28	{
29		return nsR/sR;
30	}
31	else
32	{
33		return "Infinite";
34	}
35}
36
37function BuildCodonFrequencies (obsF)
38{
39	PIStop = 1.0;
40	result = {ModelMatrixDimension,1};
41	hshift = 0;
42
43	for (h=0; h<64; h=h+1)
44	{
45		first = h$16;
46		second = h%16$4;
47		third = h%4;
48		if (_Genetic_Code[h]==10)
49		{
50			hshift = hshift+1;
51			PIStop = PIStop-obsF[first][0]*obsF[second][0]*obsF[third][0];
52			continue;
53		}
54		result[h-hshift][0]=obsF[first][0]*obsF[second][0]*obsF[third][0];
55	}
56	return result*(1.0/PIStop);
57}
58
59#include "modelParameters.mdl";
60
61HarvestFrequencies (observedFreq,filteredData,1,1,0);
62
63NICETY_LEVEL = 3;
64
65if (modelType>0)
66{
67	global R = .5;
68
69	if (modelType == 2)
70	{
71		#include "defineGamma.mdl";
72	}
73	if (modelType == 3)
74	{
75		#include "defineHM.mdl";
76	}
77}
78
79/* defines a sparse transition probabilities matrix
80 now we'll go through the matrix and assign the elements based on syn/non-syn status*/
81
82function PopulateModelMatrix (ModelMatrixName&, EFV)
83{
84	if (!ModelMatrixDimension)
85	{
86		ModelMatrixDimension = 64;
87		for (h = 0 ;h<64; h=h+1)
88		{
89			if (_Genetic_Code[h]==10)
90			{
91				ModelMatrixDimension = ModelMatrixDimension-1;
92			}
93		}
94	}
95	ModelMatrixName = {ModelMatrixDimension,ModelMatrixDimension};
96
97	hshift = 0;
98
99	if (modelType == 0)
100	{
101		for (h=0; h<64; h=h+1)
102		{
103			if (_Genetic_Code[h]==10)
104			{
105				hshift = hshift+1;
106				continue;
107			}
108			vshift = hshift;
109			for (v = h+1; v<64; v=v+1)
110			{
111				diff = v-h;
112				if (_Genetic_Code[v]==10)
113				{
114					vshift = vshift+1;
115					continue;
116				}
117			  	if ((h$4==v$4)||((diff%4==0)&&(h$16==v$16))||(diff%16==0))
118			  	{
119			  		if (h$4==v$4)
120			  		{
121			  			transition = v%4;
122			  			transition2= h%4;
123			  		}
124			  		else
125			  		{
126			  			if(diff%16==0)
127			  			{
128			  				transition = v$16;
129			  				transition2= h$16;
130			  			}
131			  			else
132			  			{
133			  				transition = v%16$4;
134			  				transition2= h%16$4;
135			  			}
136			  		}
137			  		if (_Genetic_Code[0][h]==_Genetic_Code[0][v])
138			  		{
139			  			ModelMatrixName[h-hshift][v-vshift] := synRate*EFV__[transition__][0];
140			  			ModelMatrixName[v-vshift][h-hshift] := synRate*EFV__[transition2__][0];
141				  	}
142			  		else
143			  		{
144				  		ModelMatrixName[h-hshift][v-vshift] := nonSynRate*EFV__[transition__][0];
145			  			ModelMatrixName[v-vshift][h-hshift] := nonSynRate*EFV__[transition2__][0];
146		  			}
147			  	}
148			}
149		}
150	}
151	else
152	{
153		if (modelType == 1)
154		{
155			for (h=0; h<64; h=h+1)
156				{
157					if (_Genetic_Code[h]==10)
158					{
159						hshift = hshift+1;
160						continue;
161					}
162					vshift = hshift;
163					for (v = h+1; v<64; v=v+1)
164					{
165						diff = v-h;
166						if (_Genetic_Code[v]==10)
167						{
168							vshift = vshift+1;
169							continue;
170						}
171					  	if ((h$4==v$4)||((diff%4==0)&&(h$16==v$16))||(diff%16==0))
172					  	{
173					  		if (h$4==v$4)
174					  		{
175					  			transition = v%4;
176					  			transition2= h%4;
177					  		}
178					  		else
179					  		{
180					  			if(diff%16==0)
181					  			{
182					  				transition = v$16;
183					  				transition2= h$16;
184					  			}
185					  			else
186					  			{
187					  				transition = v%16$4;
188					  				transition2= h%16$4;
189					  			}
190					  		}
191					  		if (_Genetic_Code[0][h]==_Genetic_Code[0][v])
192					  		{
193					  			ModelMatrixName[h-hshift][v-vshift] := synRate*EFV__[transition__][0];
194					  			ModelMatrixName[v-vshift][h-hshift] := synRate*EFV__[transition2__][0];
195						  	}
196					  		else
197					  		{
198						  		ModelMatrixName[h-hshift][v-vshift] := R*synRate*EFV__[transition__][0];
199					  			ModelMatrixName[v-vshift][h-hshift] := R*synRate*EFV__[transition2__][0];
200				  			}
201					  	}
202					  }
203				}
204		}
205		else
206		{
207			for (h=0; h<64; h=h+1)
208			{
209				if (_Genetic_Code[h]==10)
210				{
211					hshift = hshift+1;
212					continue;
213				}
214				vshift = hshift;
215				for (v = h+1; v<64; v=v+1)
216				{
217					diff = v-h;
218					if (_Genetic_Code[v]==10)
219					{
220						vshift = vshift+1;
221						continue;
222					}
223				  	if ((h$4==v$4)||((diff%4==0)&&(h$16==v$16))||(diff%16==0))
224				  	{
225				  		if (h$4==v$4)
226				  		{
227				  			transition = v%4;
228				  			transition2= h%4;
229				  		}
230				  		else
231				  		{
232				  			if(diff%16==0)
233				  			{
234				  				transition = v$16;
235				  				transition2= h$16;
236				  			}
237				  			else
238				  			{
239				  				transition = v%16$4;
240				  				transition2= h%16$4;
241				  			}
242				  		}
243				  		if (_Genetic_Code[0][h]==_Genetic_Code[0][v])
244				  		{
245				  			ModelMatrixName[h-hshift][v-vshift] := synRate*c*EFV__[transition__][0];
246				  			ModelMatrixName[v-vshift][h-hshift] := synRate*c*EFV__[transition2__][0];
247					  	}
248				  		else
249				  		{
250					  		ModelMatrixName[h-hshift][v-vshift] := R*synRate*c*EFV__[transition__][0];
251				  			ModelMatrixName[v-vshift][h-hshift] := R*synRate*c*EFV__[transition2__][0];
252			  			}
253				   }
254			   }
255		    }
256		}
257	}
258
259	return 0;
260}
261
262MG94 = 0;
263
264MULTIPLY_BY_FREQS = PopulateModelMatrix ("MG94", observedFreq);
265
266FREQUENCY_SENSITIVE = 1;
267
268vectorOfFrequencies = BuildCodonFrequencies (observedFreq);
269
270Model MG94model = (MG94,vectorOfFrequencies,0);
271
272if (modelType == 0)
273{
274	IS_DNDS_AVAILABLE = 1;
275}
276
277
278