1/* This file defines the transition matrix for the Muse-Gaut 94 model x an acceptance function
2
3   05/09/2005  by Sergei L. Kosakovsky Pond
4*/
5
6ModelMatrixDimension = 0;
7
8/*---------------------------------------------------------------------------------------------------------------------------------------------*/
9
10function PopulateModelMatrix (ModelMatrixName&, EFV)
11{
12	global exMean = 1;
13
14	if (!ModelMatrixDimension)
15	{
16		ModelMatrixDimension = 64;
17		for (h = 0 ;h<64; h=h+1)
18		{
19			if (_Genetic_Code[h]==10)
20			{
21				ModelMatrixDimension = ModelMatrixDimension-1;
22			}
23		}
24	}
25
26	ModelMatrixName = {ModelMatrixDimension,ModelMatrixDimension};
27
28	hshift = 0;
29
30	modelDefString = "";
31	modelDefString*16384;
32
33	if (modelType > 0)
34	{
35		catCounterAL = {};
36	}
37
38	for (h=0; h<64; h=h+1)
39	{
40		if (_Genetic_Code[h]==10)
41		{
42			hshift = hshift+1;
43			continue;
44		}
45		vshift = hshift;
46		for (v = h+1; v<64; v=v+1)
47		{
48			diff = v-h;
49			if (_Genetic_Code[v]==10)
50			{
51				vshift = vshift+1;
52				continue;
53			}
54			nucPosInCodon = 2;
55			if ((h$4==v$4)||((diff%4==0)&&(h$16==v$16))||(diff%16==0))
56			{
57				if (h$4==v$4)
58				{
59					transition = v%4;
60					transition2= h%4;
61				}
62				else
63				{
64					if(diff%16==0)
65					{
66						transition = v$16;
67						transition2= h$16;
68						nucPosInCodon = 0;
69					}
70					else
71					{
72						transition = v%16$4;
73						transition2= h%16$4;
74						nucPosInCodon = 1;
75					}
76				}
77				hs = Format(h-hshift,0,0);
78				vs = Format(v-vshift,0,0);
79				ts = Format(transition,0,0);
80				ts2= Format(transition2,0,0);
81				ps = Format(nucPosInCodon,0,0);
82				aa1 = _Genetic_Code[0][h];
83				aa2 = _Genetic_Code[0][v];
84				if (aa1==aa2)
85				{
86					modelDefString*("ModelMatrixName["+hs+"]["+vs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*EFV__["+ts+"]["+ps+"];\n"+
87													 "ModelMatrixName["+vs+"]["+hs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*EFV__["+ts2+"]["+ps+"];\n");
88				}
89				else
90				{
91					bt = aaRateMultipliers[aa1][aa2];
92					if (modelType > 0)
93					{
94						modelDefString*("ModelMatrixName["+hs+"]["+vs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*"+bt+"^(exMean*c)*EFV__["+ts+"]["+ps+"];\n"+
95														 "ModelMatrixName["+vs+"]["+hs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*"+bt+"^(exMean*c)*EFV__["+ts2+"]["+ps+"];\n");
96					}
97					else
98					{
99						bt = aaRateMultipliers[aa1][aa2];
100						modelDefString*("ModelMatrixName["+hs+"]["+vs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*"+bt+"^exMean*EFV__["+ts+"]["+ps+"];\n"+
101														 "ModelMatrixName["+vs+"]["+hs+"] := "+_nucBiasTerms[transition][transition2]+"synRate*"+bt+"^exMean*EFV__["+ts2+"]["+ps+"];\n");
102					}
103				}
104			}
105	    }
106    }
107	modelDefString*0;
108	ExecuteCommands (modelDefString);
109	return 0;
110}
111
112
113/*---------------------------------------------------------------------------------------------------------------------------------------------*/
114
115function BuildCodonFrequencies (obsF)
116{
117	PIStop = 1.0;
118	result = {ModelMatrixDimension,1};
119	hshift = 0;
120
121	for (h=0; h<64; h=h+1)
122	{
123		first = h$16;
124		second = h%16$4;
125		third = h%4;
126		if (_Genetic_Code[h]==10)
127		{
128			hshift = hshift+1;
129			PIStop = PIStop-obsF[first][0]*obsF[second][1]*obsF[third][2];
130			continue;
131		}
132		result[h-hshift][0]=obsF[first][0]*obsF[second][1]*obsF[third][2];
133	}
134	return result*(1.0/PIStop);
135}
136
137/*---------------------------------------------------------------------------------------------------------------------------------------------*/
138
139categoriesUsed = 0;
140
141#include "MGwEX.ibf";
142#include "modelParameters2.mdl";
143
144global AC = 1;
145global AT = 1;
146global CG = 1;
147global CT = 1;
148global GT = 1;
149
150sharedFlag = 1;
151if (modelType >0)
152{
153	if (modelType == 1)
154	{
155		categoriesUsed = 1;
156		#include "defineGamma.mdl";
157	}
158	if (modelType == 2)
159	{
160		categoriesUsed = 1;
161		#include "defineHM.mdl";
162	}
163}
164
165MGCustomRateBiasTerms = {{"AC*","","AT*","CG*","CT*","GT*"}};
166
167done = 0;
168while (!done)
169{
170	fprintf (stdout,"\nPlease enter a 6 character model designation (e.g:010010 defines HKY85):");
171	fscanf  (stdin,"String", modelDesc);
172	if (Abs(modelDesc)==6)
173	{
174		done = 1;
175	}
176}
177
178
179paramCount	  = 0;
180_nucBiasTerms = {4,4};
181_nucBiasTerms[0][0] = "";
182
183
184if (modelDesc[0]==modelDesc[1])
185{
186	MGCustomRateBiasTerms[0] = MGCustomRateBiasTerms[1];
187}
188
189_nucBiasTerms[1][0] = MGCustomRateBiasTerms[0];
190_nucBiasTerms[0][1] = MGCustomRateBiasTerms[0];
191_nucBiasTerms[2][0] = MGCustomRateBiasTerms[1];
192_nucBiasTerms[0][2] = MGCustomRateBiasTerms[1];
193
194h = 0;
195v = 3;
196
197for (customLoopCounter2=2; customLoopCounter2<6; customLoopCounter2=customLoopCounter2+1)
198{
199	for (customLoopCounter=0; customLoopCounter<customLoopCounter2; customLoopCounter=customLoopCounter+1)
200	{
201		if (modelDesc[customLoopCounter]==modelDesc[customLoopCounter2])
202		{
203			_nucBiasTerms[h][v] = MGCustomRateBiasTerms[customLoopCounter];
204			_nucBiasTerms[v][h] = MGCustomRateBiasTerms[customLoopCounter];
205			break;
206		}
207	}
208	if (customLoopCounter == customLoopCounter2)
209	{
210		_nucBiasTerms[h][v] = MGCustomRateBiasTerms[customLoopCounter2];
211		_nucBiasTerms[v][h] = MGCustomRateBiasTerms[customLoopCounter2];
212	}
213
214	v = v+1;
215	if (v==4)
216	{
217		h=h+1;
218		v=h+1;
219	}
220}
221
222
223HarvestFrequencies (observedFreq,filteredData,3,1,1);
224
225
226NICETY_LEVEL = 3;
227
228MG94custom = 0;
229
230MULTIPLY_BY_FREQS = PopulateModelMatrix ("MG94custom", observedFreq);
231
232FREQUENCY_SENSITIVE = 1;
233
234vectorOfFrequencies = BuildCodonFrequencies (observedFreq);
235
236Model MG94customModel = (MG94custom,vectorOfFrequencies,0);
237
238USE_POSITION_SPECIFIC_FREQS = 1;
239