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