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