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