1skipCodeSelectionStep = 0; 2 3 4LoadFunctionLibrary("chooseGeneticCode"); 5LoadFunctionLibrary("GrabBag"); 6LoadFunctionLibrary("dSdNTreeTools"); 7LoadFunctionLibrary("CF3x4"); 8LoadFunctionLibrary("BranchSiteTemplate"); 9LoadFunctionLibrary("DescriptiveStatistics"); 10LoadFunctionLibrary("TreeTools"); 11 12 13DataSet ds = ReadDataFile(PROMPT_FOR_FILE); 14DataSetFilter dsf = CreateFilter(ds,3,"","",GeneticCodeExclusions); 15 16HarvestFrequencies (nuc3, dsf, 3, 1, 1); 17nucCF = CF3x4 (nuc3, GeneticCodeExclusions); 18 19rate_class_count = prompt_for_a_value ("How many rate classes", 3, 2, 8, 1); 20fprintf (stdout, "\nUsing ", rate_class_count, " rate classes\n"); 21 22PopulateModelMatrix ("MGMatrixLocal", nucCF, "syn", "", "nonsyn"); 23codon3x4 = BuildCodonFrequencies (nucCF); 24Model MGL = (MGMatrixLocal, codon3x4, 0); 25 26_DO_TREE_REBALANCE_ = 0; 27LoadFunctionLibrary ("queryTree"); 28 29SetDialogPrompt ("Save analysis results to"); 30fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN,"Model\tLogL\tNP\tBIC\tTree String"); 31csvFilePath = LAST_FILE_PATH; 32 33 34// ---- PHASE 0 : LOCAL (no omega variation model) FIT. ----- 35 36fprintf (stdout, "[BS-REL PHASE 0] Fitting the local MG94 (no site-to-site variation) to obtain initial parameter estimates\n"); 37 38 39LikelihoodFunction base_LF = (dsf, givenTree); 40 41/* 42ExecuteAFile (csvFilePath + ".mglocal.fit"); 43res_base = {2,2}; 44*/ 45 46Optimize (res_base,base_LF); 47writeTheLF (".mglocal.fit", "base_LF"); 48baseParameters = 9; 49localLL = res_base[1][0]; 50localParams = res_base[1][1] + baseParameters; 51 52totalBranchCount = BranchCount(givenTree) + TipCount (givenTree); 53 54pValueByBranch = {totalBranchCount,10}; 55bNames = BranchName (givenTree, -1); 56 57for (k = 0; k < totalBranchCount; k = k+1) { 58 srate = Eval ("givenTree." + bNames[k] + ".syn"); 59 nsrate = Eval ("givenTree." + bNames[k] + ".nonsyn"); 60 if (srate > 0) { 61 pValueByBranch [k][0] = Min (10, nsrate/srate); 62 } else { 63 pValueByBranch [k][0] = 10; 64 } 65} 66 67sample_size = dsf.sites * dsf.species; 68omegaStats = GatherDescriptiveStats (pValueByBranch[-1][0]); 69fprintf (stdout, "Local omega model (no rate variation) "); 70 71 72localBIC = BIC (localLL, localParams, sample_size); 73PrintDescriptiveStats ("Branch omega values", omegaStats); 74 75for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 76 PopulateModelMatrix ("MGMatrix"+rcc, nucCF, "t", "omega"+rcc, ""); 77 ExecuteCommands ("global omegaG" + rcc + " = " + (0.3*rcc - 0.2)); 78 PopulateModelMatrix ("MGMatrixG"+rcc, nucCF, "t", "omegaG"+rcc, ""); 79} 80 81// ---- PHASE 1 : FULL GLOBAL MODEL FIT ---- // 82 83generate_bs_rel_gdd_freqs (rate_class_count, "global_freqsG", "model_mixingG", "PauxG", "MGMatrixG",1); 84generate_bs_rel_gdd_freqs (rate_class_count, "global_freqs", "model_mixing", "Paux", "MGMatrix",0); 85 86sortedLocalDNDS = (pValueByBranch[-1][0])%0; 87 88for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 89 ExecuteCommands ("omegaG" + rcc + "=" + sortedLocalDNDS[Max(0,totalBranchCount*rcc$rate_class_count-1)]); 90} 91 92 93ExecuteCommands ("Model MGG =(\"`model_mixingG`\",codon3x4,EXPLICIT_FORM_MATRIX_EXPONENTIAL);"); 94Tree mixtureTreeG = treeString; 95 96ExecuteCommands ("Model MGLocal=(\"`model_mixing`\",codon3x4,EXPLICIT_FORM_MATRIX_EXPONENTIAL);"); 97Tree mixtureTree = treeString; 98 99ReplicateConstraint ("this1.?.t:=this2.?.syn",mixtureTreeG,givenTree); 100ReplicateConstraint ("this1.?.t:=this2.?.syn",mixtureTree,givenTree); 101ClearConstraints (mixtureTree); 102ClearConstraints (mixtureTreeG); 103 104 105ASSUME_REVERSIBLE_MODELS = 1; 106OPTIMIZATION_METHOD = 4; 107 108BIC_scores = {}; 109BranchLengthEstimates = {}; 110 111mixTreeAVL = mixtureTreeG ^ 0; 112 113LikelihoodFunction three_LFG = (dsf,mixtureTreeG); 114fprintf (stdout, "[BS-REL PHASE 1] Fitting a GLOBAL branch-site matrix mixture\n"); 115 116 117blG = getBranchLengthExpression (1); 118removeT = "/(" + (blG ^{{"\\*t",""}}) + ")"; 119 120fprintf (stdout, base_LF); 121 122blLocal = {}; 123for (k = 0; k < totalBranchCount; k += 1) { 124 blLocal [bNames[k]] = Eval("BranchLength(givenTree,\""+bNames[k]+"\")")*3; 125 ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t := " + blLocal [bNames[k]] + removeT ); 126} 127 128 129USE_LAST_RESULTS = 1; 130SHORT_MPI_RETURN = 1; 131VERBOSITY_LEVEL = 0; 132 133OPTIMIZATION_PRECISION = 0.1; 134 135Optimize (res_three_LF_global,three_LFG); 136 137fprintf (stdout, "\n\nApproximation phase 1 (use MG-local branch lengths): ", res_three_LF_global[1][0], "\n"); 138fprintf (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,getBranchLengths ("mixtureTreeG", 1))); 139 140OPTIMIZATION_PRECISION = 0.01; 141 142for (k = 0; k < totalBranchCount; k += 1) { 143 ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t = mixtureTreeG." + bNames[k] + ".t"); 144 bl = blLocal [bNames[k]] ; 145 ExecuteCommands ("FindRoot (z,`blG`-"+bl+",t,0,10000);"); 146 ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t :< " + 10*z); 147 148} 149 150Optimize (res_three_LF_global,three_LFG); 151fprintf (stdout, "\n\nApproximation phase 2 (maximum branch lengths are limited to 10x those of the null model): ", res_three_LF_global[1][0], "\n"); 152fprintf (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,getBranchLengths ("mixtureTreeG", 1))); 153 154OPTIMIZATION_PRECISION = 0.001; 155 156for (k = 0; k < totalBranchCount; k += 1) { 157 ExecuteCommands ("mixtureTreeG." + bNames[k] + ".t :< 10000"); 158} 159 160Optimize (res_three_LF_global,three_LFG); 161 162//ExecuteAFile (csvFilePath + ".relglobal.fit"); 163 164 165blStashByName = {}; 166for (k = 0; k < totalBranchCount; k += 1) { 167 bl = blLocal [bNames[k]] ; 168 ExecuteCommands ("FindRoot (z,`blG`-"+bl+",t,0,10000);"); 169 blStashByName[bNames[k]] = z; 170 //fprintf (stdout, "\n", bNames[k], " : ", z, " (", bl/3, ")\n"); 171} 172 173writeTheLF (".relglobal.fit", "three_LFG"); 174 175 176fprintf (stdout, "Global model fit:"); 177 178BIC_scores ["global model"] = BIC (res_three_LF_global[1][0], res_three_LF_global[1][1]+baseParameters, sample_size); 179BranchLengthEstimates ["global model"] = getBranchLengths ("mixtureTreeG", 1); 180 181fprintf (csvFilePath, "\nGlobal\t", res_three_LF_global[1][0],"\t", res_three_LF_global[1][1]+baseParameters, "\t", BIC_scores ["global model"], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["global model"])); 182 183fprintf (stdout, "\nInferred this global omega distribution: "); 184reportOmegaDistro (""); 185fprintf (stdout,"Tree branch lengths:\n",PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["global model"])); 186globalState = saveLF ("three_LFG"); 187 188 189for (k = 0; k < totalBranchCount; k += 1) { 190 constrainABranch (bNames[k], 0); 191} 192 193// ---- PHASE 2 : LOCAL MODEL FITS ---- // 194 195fprintf (stdout, "\n[BS-REL PHASE 2] Fitting a GLOBAL branch-site matrix mixture with a SINGLE unconstrained branch\n"); 196 197OPTIMIZATION_METHOD = 0; 198 199LikelihoodFunction three_LF = (dsf,mixtureTree); 200 201branchValues = {}; 202 203if (MPI_NODE_COUNT > 1) { 204 MPI_NODE_STATE = {MPI_NODE_COUNT-1,1}; 205 MPI_NODE_STATE[0] = ""; 206} 207 208 209for (k = 0; k < totalBranchCount; k+=1) { 210 thisBranchName = bNames[k]; 211 fprintf (stdout, "\n[BS-REL PHASE 2. Branch '", thisBranchName, "']\n"); 212 globalState["restoreLF"][""]; 213 if (k) { 214 if (MPI_NODE_COUNT > 1) { 215 for (nodeID = 0; nodeID < totalBranchCount; nodeID += 1) { 216 constrainABranch (bNames[nodeID],0); 217 } 218 } else { 219 constrainABranch (bNames[k-1],0); 220 } 221 } 222 unConstrainABranch (thisBranchName); 223 224 if (MPI_NODE_COUNT > 1) { 225 for (nodeID = 0; nodeID < MPI_NODE_COUNT-1; nodeID += 1) { 226 if (Abs(MPI_NODE_STATE[nodeID]) == 0) { 227 MPISend (nodeID+1, three_LF); 228 fprintf (stdout, "\n[SENT TO NODE ", nodeID, "]\n"); 229 MPI_NODE_STATE [nodeID] = thisBranchName; 230 break; 231 } 232 } 233 if (nodeID == MPI_NODE_COUNT-1) { 234 processABranch (thisBranchName, 1); 235 } 236 } else { 237 OPTIMIZATION_PRECISION = 0.05; 238 maxBL = blStashByName[thisBranchName]*10; 239 t = maxBL; 240 fprintf (stdout, thisBranchName, ": maximum t value: ", maxBL, ". Maximum achieved branch length: ", Eval (blG), "\n"); 241 constrainABranch (thisBranchName,maxBL); 242 Optimize (localBranchRes, three_LF); 243 fprintf (stdout, "Approximation phase 1 (constrained branch length): ", localBranchRes[1][0], "\n"); 244 constrainABranch (thisBranchName,-1); 245 OPTIMIZATION_PRECISION = 0.001; 246 Optimize (localBranchRes, three_LF); 247 processABranch (thisBranchName,0); 248 } 249} 250 251 252leftOver = 0; 253for (nodeID = 0; nodeID < MPI_NODE_COUNT-1; nodeID += 1) { 254 leftOver += Abs(MPI_NODE_STATE[nodeID])>0; 255} 256 257//fprintf (stdout, MPI_NODE_STATE, "\n", leftOver, "\n"); 258 259 260for (nodeID = 0; nodeID < leftOver; nodeID += 1) { 261 processABranch ("", 0); 262} 263 264branchValues ["restoreLF"][""]; 265 266fprintf (stdout, "\n[BS-REL PHASE 3] Fitting a LOCAL branch-site matrix mixture model\n"); 267Optimize (res_local, three_LF); 268 269writeTheLF (".local.fit", "three_LF"); 270 271fprintf (stdout, "Local model fit:"); 272BIC_scores ["local model"] = BIC (res_local[1][0], res_local[1][1]+baseParameters, sample_size); 273BranchLengthEstimates ["local model"] = getBranchLengths ("mixtureTree", 0); 274fprintf (stdout, "\nTree: ", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["local model"]), "\n"); 275fprintf (csvFilePath, "\nLocal\t", res_local[1][0],"\t", res_local[1][1]+baseParameters, "\t", BIC_scores ["local model"], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates ["local model"]), CLOSE_FILE); 276 277min_BIC = 1e100; 278 279function compute_min_BIC (key, value) { 280 min_BIC = Min (min_BIC, value); 281 return 0; 282} 283 284function akaike_weights_pass1 (key, value) { 285 akaike_scores [key] = Exp (0.5*(min_BIC-value)); 286 return 0; 287} 288 289 290function akaike_weights_pass2 (key, value) { 291 BIC_scores [key] = value/norm_factor; 292 return 0; 293} 294 295fprintf (stdout, "\n[BS-REL PHASE 4] Generating a model averaged branch lengths\n"); 296akaike_scores = {}; 297BIC_scores ["compute_min_BIC"][""]; 298BIC_scores ["akaike_weights_pass1"][""]; 299norm_factor = +akaike_scores; 300akaike_scores ["akaike_weights_pass2"][""]; 301 302reweighted_branch_lengths = {}; 303 304model_names = Rows(BranchLengthEstimates); 305for (k = 0; k < Abs (BranchLengthEstimates); k+=1) { 306 key = model_names [k]; 307 bl_avl = BranchLengthEstimates[key]; 308 bnames = Rows (bl_avl); 309 310 for (b = 0; b < Abs (bl_avl); b += 1) { 311 reweighted_branch_lengths[bnames[b]] += BIC_scores[key]*bl_avl[bnames[b]]; 312 } 313} 314 315fprintf (stdout, PostOrderAVL2StringDistances(mixTreeAVL,reweighted_branch_lengths), "\n"); 316 317//------------------------------------------------------------------------------------------------------------------------ 318function processABranch (thisBranchName, doSend) { 319 if (MPI_NODE_COUNT > 1) { 320 MPIReceive (-1,fromNode,resStr); 321 prevBranch = MPI_NODE_STATE[fromNode-1]; 322 fprintf (stdout, "\n[RECEIVED " + prevBranch +" FROM NODE ", (fromNode-1), "]\n"); 323 if (doSend) { 324 MPISend (fromNode, three_LF); 325 MPI_NODE_STATE [fromNode-1] = thisBranchName; 326 fprintf (stdout, "\n[SENT " + thisBranchName +" TO NODE ", (fromNode-1), "]\n"); 327 } else { 328 MPI_NODE_STATE [fromNode-1] = ""; 329 } 330 thisBranchName = prevBranch; 331 ExecuteCommands (resStr); 332 three_LF_MLE_VALUES ["restoreLF"][""]; 333 localBranchRes = three_LF_MLES; 334 } 335 writeTheLF (".`thisBranchName`.fit", "three_LF"); 336 fprintf (stdout, "\nModel fit:"); 337 BIC_scores [thisBranchName] = BIC (localBranchRes[1][0], localBranchRes[1][1]+baseParameters, sample_size); 338 BranchLengthEstimates [thisBranchName] = getBranchLengths ("mixtureTree", 0); 339 fprintf (csvFilePath, "\nGlobal+",thisBranchName,"\t", localBranchRes[1][0],"\t", localBranchRes[1][1]+baseParameters, "\t", BIC_scores [thisBranchName], "\t", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates [thisBranchName])); 340 fprintf (stdout, "\nLocal branch omega distribution: "); 341 reportOmegaDistro(thisBranchName); 342 fprintf (stdout, "Global omega distribution on the rest of the branches: "); 343 reportOmegaDistro (""); 344 stashBranchValues (thisBranchName, "branchValues"); 345 pv = 1-CChi2 (2*(localBranchRes[1][0]-res_three_LF_global[1][0]),5); 346 fprintf (stdout, "\nLRT p-value for branch deviation from the global pattern = ", pv, "\n"); 347 fprintf (stdout,"\n", PostOrderAVL2StringDistances(mixTreeAVL,BranchLengthEstimates [thisBranchName])); 348 return 0; 349 350} 351 352//------------------------------------------------------------------------------------------------------------------------ 353function writeTheLF (fileNameExt,lfID) { 354 lfOut = csvFilePath + fileNameExt;//".local.fit"; 355 LIKELIHOOD_FUNCTION_OUTPUT = 7; 356 ExecuteCommands ("fprintf (lfOut, CLEAR_FILE, `lfID`);"); 357 fprintf (stdout, "[WROTE MODEL FIT TO ", lfOut, "]\n"); 358 LIKELIHOOD_FUNCTION_OUTPUT = 2; 359 return 0; 360} 361 362//------------------------------------------------------------------------------------------------------------------------ 363function constrainABranch (branch_name, z) { 364 if (z > 0) { 365 ExecuteCommands ("mixtureTree." + branch_name + ".t :< " + z); 366 return 0; 367 } else { 368 if (z < 0) { 369 ExecuteCommands ("mixtureTree." + branch_name + ".t :< 10000"); 370 return 0; 371 } 372 } 373 ExecuteCommands ("mixtureTree." + branch_name + ".t = mixtureTreeG." + branch_name + ".t"); 374 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 375 ExecuteCommands ("mixtureTree." + branch_name + ".omega" + rcc + ":= omegaG" + rcc); 376 } 377 378 for (rcc = 1; rcc < rate_class_count; rcc += 1) { 379 ExecuteCommands ("mixtureTree." + branch_name + ".Paux" + rcc + ":= PauxG" + rcc); 380 } 381 382 return 0; 383} 384 385//------------------------------------------------------------------------------------------------------------------------ 386function stashBranchValues (branch_name, storage&) { 387 388 locals = {"0" : "t"}; 389 390 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 391 locals + ("omega" + rcc); 392 if (rcc < rate_class_count) { 393 locals + ("Paux" + rcc); 394 } 395 } 396 for (_varID = 0; _varID < Abs (locals); _varID += 1) { 397 storage ["mixtureTree." + branch_name + "." + locals[_varID]] = Eval ("mixtureTree." + branch_name + "." + locals[_varID]); 398 } 399 400 return 0; 401} 402 403//------------------------------------------------------------------------------------------------------------------------ 404function unConstrainABranch (branch_name) { 405 406 ExecuteCommands ("mixtureTree." + branch_name + ".t = mixtureTreeG." + branch_name + ".t"); 407 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 408 ExecuteCommands ("mixtureTree." + branch_name + ".omega" + rcc + "= omegaG" + rcc); 409 } 410 411 for (rcc = 1; rcc < rate_class_count; rcc += 1) { 412 ExecuteCommands ("mixtureTree." + branch_name + ".Paux" + rcc + "= PauxG" + rcc); 413 } 414 415 return 0; 416} 417 418//------------------------------------------------------------------------------------------------------------------------ 419function BIC (logL, params, sample_size) { 420 BICvalue = -2*logL + 2*params*(sample_size-params-1); 421 fprintf (stdout, " log(L) = ", logL, " with ", params, " parameters, yielding BIC = ", BICvalue, " assuming the sample size of ", sample_size, "."); 422 return BICvalue; 423} 424 425//------------------------------------------------------------------------------------------------------------------------ 426function reportOmegaDistro (branchName) { 427 if (Abs(branchName)==0) { 428 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 429 fprintf (stdout, "\n\t omega = ", Eval("omegaG"+rcc), ", p = ", Eval(global_freqsG[rcc-1])); 430 } 431 } else { 432 for (rcc = 1; rcc < rate_class_count; rcc += 1) { 433 ExecuteCommands ("Paux" + rcc + "=mixtureTree.`branchName`.Paux" +rcc); 434 } 435 436 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 437 fprintf (stdout, "\n\t omega = ", Eval("mixtureTree.`branchName`.omega"+rcc), ", p = ", Eval(global_freqs[rcc-1])); 438 } 439 } 440 fprintf (stdout, "\n"); 441 return 0; 442} 443 444//------------------------------------------------------------------------------------------------------------------------ 445 446function getBranchLengthExpression (modelType) { 447 bl_expression = ""; bl_expression * 128; 448 449 if (modelType == 0) { 450 modelPrefix = "MGMatrix"; 451 freqMatrix = global_freqs; 452 } else { 453 modelPrefix = "MGMatrixG"; 454 freqMatrix = global_freqsG; 455 } 456 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 457 if (rcc > 1) { 458 bl_expression * "+"; 459 } 460 ExecuteCommands ("Model _temp = (`modelPrefix`" + rcc + ", codon3x4, 0);"); 461 GetString (modelBL, _temp, -1); 462 bl_expression * (freqMatrix[rcc-1]+"*(`modelBL`)"); 463 } 464 bl_expression * 0; 465 return bl_expression; 466} 467 468//------------------------------------------------------------------------------------------------------------------------ 469 470function getBranchLengths (treeID, modelType) { 471 blByName = {}; 472 bnames_list = Eval ("BranchName (`treeID`,-1)"); 473 474 bl_expression = getBranchLengthExpression(modelType); 475 476 if (modelType == 0) { 477 locals = {"0" : "t"}; 478 479 for (rcc = 1; rcc <= rate_class_count; rcc += 1) { 480 locals + ("omega" + rcc); 481 if (rcc < rate_class_count) { 482 locals + ("Paux" + rcc); 483 } 484 } 485 486 for (_bID = 0; _bID < Columns (bnames_list) - 1; _bID += 1) { 487 branch_name = bnames_list[_bID]; 488 489 for (_varID = 0; _varID < Abs (locals); _varID += 1) { 490 ExecuteCommands (locals[_varID] + " = Eval (\"`treeID`.`branch_name`." + locals[_varID] + "\")"); 491 } 492 blByName [branch_name] = Eval (bl_expression)/3; 493 } 494 } else { 495 for (_bID = 0; _bID < Columns (bnames_list) - 1; _bID += 1) { 496 branch_name = bnames_list[_bID]; 497 t = Eval ("`treeID`.`branch_name`.t"); 498 blByName [branch_name] = Eval (bl_expression)/3; 499 } 500 } 501 502 return blByName; 503} 504 505 506//------------------------------------------------------------------------------------------------------------------------ 507function saveLF (ID) 508{ 509 ExecuteCommands ("GetString(_lfInfo,"+ID+",-1)"); 510 _stashLF = {}; 511 for (_k = 0; _k < Columns (_lfInfo["Global Independent"]); _k+=1) 512 { 513 _stashLF [(_lfInfo["Global Independent"])[_k]] = Eval ((_lfInfo["Global Independent"])[_k]); 514 } 515 for (_k = 0; _k < Columns (_lfInfo["Local Independent"]); _k+=1) 516 { 517 _stashLF [(_lfInfo["Local Independent"])[_k]] = Eval ((_lfInfo["Local Independent"])[_k]); 518 } 519 520 return _stashLF; 521} 522 523//------------------------------------------------------------------------------------------------------------------------ 524function restoreLF (key, value) 525{ 526 ExecuteCommands (key + " = " + value); 527 return 0; 528} 529 530 531//------------------------------------------------------------------------------------------------------------------------ 532function generate_bs_rel_gdd_freqs (numberOfRates, freqs&, model_mixing&, probPrefix, modelPrefix, is_global) 533{ 534 535 freqs = {numberOfRates,1}; 536 model_mixing = ""; model_mixing * 128; 537 538 if (numberOfRates == 1) { 539 freqs[0] = "1"; 540 } 541 else { 542 global_prefix = ""; 543 544 if (is_global) 545 { global_prefix = "global "; } 546 547 for (mi=1; mi<numberOfRates; mi += 1) { 548 ExecuteCommands (global_prefix+probPrefix+mi+":<1;"+probPrefix+mi+" = 1/" + (numberOfRates-mi+1)); 549 ExecuteCommands (global_prefix+probPrefix+mi+":>0;"); 550 } 551 552 freqs[0] = ""+probPrefix+"1"; 553 for (mi=1; mi<numberOfRates-1; mi+=1) { 554 freqs[mi] = ""; 555 for (mi2=1;mi2<=mi;mi2+=1) { 556 freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi2+")*"; 557 } 558 freqs[mi] = freqs[mi] + probPrefix+(mi+1); 559 } 560 561 freqs[mi] = ""; 562 for (mi2=1;mi2<mi;mi2+=1) 563 { 564 freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi2+")*"; 565 } 566 freqs[mi] = freqs[mi] + "(1-"+probPrefix+mi+")"; 567 } 568 569 model_mixing * ("Exp(`modelPrefix`1)*"+freqs[0]); 570 for (mi = 1; mi < numberOfRates; mi=mi+1) { 571 model_mixing * ("+Exp(`modelPrefix`"+(mi+1)+")*" + freqs[mi]); 572 } 573 model_mixing * 0; 574 return 0; 575} 576