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