1RequireVersion ("2.00.20100101"); 2 3ChoiceList (runType,"Run Type",1,SKIP_NONE, 4 "New run","Start a new run", 5 "Continued run","Start with an initial point output by a Gateux derivative approximator." 6 ); 7 8if (runType < 0) 9{ 10 return 0; 11} 12 13 14randomizeInitValues = 0; 15ModelMatrixDimension = 0; 16 17 18/*---------------------------------------------------------------------------------------------------------------------------------------------------*/ 19 20function ReportDistributionString (rc) 21{ 22 distroString = ""; 23 distroString * 1024; 24 25 reportMx = {rc + (bivariateFitHasMultipleCladeRates>1)*(bivariateFitHasMultipleCladeRates-1),4}; 26 27 for (mi=0; mi<rc; mi=mi+1) 28 { 29 reportMx[mi][0] = Eval ("S_"+mi+"/c_scale"); 30 reportMx[mi][1] = Eval ("NS_"+mi+"/c_scale"); 31 reportMx[mi][2] = reportMx[mi][1]/reportMx[mi][0]; 32 reportMx[mi][3] = Eval (freqStrMx[mi]); 33 34 distroString * ("Class "+(mi+1)+"\n\tdS = "+Format(reportMx[mi][0],10,3) 35 +"\n\tdN = "+Format(reportMx[mi][1],10,3) 36 +"\n\tdN/dS = "+Format(reportMx[mi][2],10,3) 37 +"\n\tProb = "+Format(reportMx[mi][3],10,3)); 38 39 if (bivariateFitHasMultipleCladeRates > 1) 40 { 41 distroString * "\n\tBranch group-pecific dN/dS"; 42 43 for (mi2 = 1; mi2 < bivariateFitHasMultipleCladeRates; mi2 += 1) 44 { 45 distroString * ( "\n\tGroup " + Format (mi2,2,0) + " dN = " + Format (reportMx[mi][1]*Eval("clade_"+mi2+"_NS_"+mi),10,3)); 46 distroString * ( "\n\tGroup " + Format (mi2,2,0) + " dN/dS = " + Format (reportMx[mi][2]*Eval("clade_"+mi2+"_NS_"+mi),10,3)); 47 } 48 } 49 50 distroString * "\n"; 51 52 } 53 54 distroString * 0; 55 return distroString; 56} 57 58 59 60totalCodonCount = 0; 61totalCharCount = 0; 62 63/*---------------------------------------------------------------------------------------------------------------------------------------------------*/ 64if (runType == 0) 65{ 66 LF_NEXUS_EXPORT_EXTRA = "bivariateFitHasMultipleCladeRates = " + bivariateFitHasMultipleCladeRates + ";"; 67 68 ChoiceList (branchLengths,"Branch Lengths",1,SKIP_NONE, 69 "Codon Model","Jointly optimize rate parameters and branch lengths (slow and thorough)", 70 "Nucleotide Model","Estimate branch lengths once, using an appropriate nucleotide model (quick and dirty)." 71 ); 72 73 if (branchLengths<0) 74 { 75 return; 76 } 77 78 fileCount = 0; 79 while (fileCount < 1) 80 { 81 fprintf (stdout, "How many datafiles are to be analyzed (>=1):?"); 82 fscanf (stdin, "Number", fileCount); 83 fileCount = fileCount $ 1; 84 } 85 86 #include "TemplateModels/chooseGeneticCode.def"; 87 88 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 89 { 90 fprintf (stdout, "\nReading input file ", fileID, "/", fileCount, "\n"); 91 SetDialogPrompt ("Please specify codon data #" + fileID + ":"); 92 93 ExecuteCommands ( 94 "DataSet ds_" +fileID + " = ReadDataFile (PROMPT_FOR_FILE);" + 95 "fprintf (stdout,\"\\n______________READ THE FOLLOWING DATA______________\\n\",ds_" + fileID + ");"+ 96 "DataSetFilter filteredData_" + fileID + " = CreateFilter (ds_"+fileID +",3,\"\",\"\",GeneticCodeExclusions);"); 97 #include "queryTree.bf"; 98 ExecuteCommands ("treeString_" + fileID + " = treeString;totalCodonCount=totalCodonCount+filteredData_"+ 99 fileID + ".sites;totalCharCount=totalCharCount+filteredData_" + fileID + ".sites*filteredData_" + fileID + ".species;"); 100 ExecuteCommands ("treeStringNoModels_" + fileID + "= treeString_" + fileID + " ^{{\"\\\\{[^\\\\}]+\\\\}\",\"\"}}"); 101 102 } 103 104 observedFreq = {4,3}; 105 observedFreqSingle = {4,1}; 106 107 108 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 109 { 110 ExecuteCommands ("HarvestFrequencies (tp, filteredData_"+fileID+",3,1,1);HarvestFrequencies (ts, filteredData_"+fileID+",1,1,1);cfs = filteredData_"+fileID+".sites;"); 111 observedFreq = observedFreq + tp*(cfs/totalCodonCount); 112 observedFreqSingle = observedFreqSingle + ts*(cfs/totalCodonCount); 113 } 114 115 ExecuteAFile ("TemplateModels/CF3x4.bf"); 116 observedFreq = CF3x4 (observedFreq, GeneticCodeExclusions); 117 118 done = 0; 119 while (!done) 120 { 121 fprintf (stdout,"\nPlease enter a 6 character model designation (e.g:010010 defines HKY85):"); 122 fscanf (stdin,"String", modelDesc); 123 if (Abs(modelDesc)==6) 124 { 125 done = 1; 126 } 127 } 128 129 ModelTitle = "MG94x"+modelDesc[0]; 130 131 rateBiasTerms = {{"AC","1","AT","CG","CT","GT"}}; 132 paramCount = 0; 133 134 modelConstraintString = ""; 135 136 for (customLoopCounter2=1; customLoopCounter2<6; customLoopCounter2=customLoopCounter2+1) 137 { 138 for (customLoopCounter=0; customLoopCounter<customLoopCounter2; customLoopCounter=customLoopCounter+1) 139 { 140 if (modelDesc[customLoopCounter2]==modelDesc[customLoopCounter]) 141 { 142 ModelTitle = ModelTitle+modelDesc[customLoopCounter2]; 143 if (rateBiasTerms[customLoopCounter2] == "1") 144 { 145 modelConstraintString = modelConstraintString + rateBiasTerms[customLoopCounter]+":="+rateBiasTerms[customLoopCounter2]+";"; 146 } 147 else 148 { 149 modelConstraintString = modelConstraintString + rateBiasTerms[customLoopCounter2]+":="+rateBiasTerms[customLoopCounter]+";"; 150 } 151 break; 152 } 153 } 154 if (customLoopCounter==customLoopCounter2) 155 { 156 ModelTitle = ModelTitle+modelDesc[customLoopCounter2]; 157 } 158 } 159 160 ExecuteAFile ("2RatesAnalyses/MG94xREVxBivariate_Multirate.mdl"); 161 162 if (Abs(modelConstraintString)) 163 { 164 ExecuteCommands (modelConstraintString); 165 } 166 167 ChoiceList (randomizeInitValues, "Initial Value Options",1,SKIP_NONE, 168 "Default", "Use default inital values for rate distribution parameters.", 169 "Randomized", "Select initial values for rate distribution parameters at random.", 170 "User Values", "Set initial values for rate distribution parameters and fix them for the optimization process"); 171 172 173 if (randomizeInitValues < 0) 174 { 175 return; 176 } 177 178 resp = 0; 179 180 while (resp<1) 181 { 182 fprintf (stdout,"Number of rate classes:"); 183 fscanf (stdin, "Number", resp); 184 } 185 186 if (randomizeInitValues < 2) 187 { 188 ChoiceList (stratDNDS, "dN/dS Stratification",1,SKIP_NONE, 189 "Unconstrained", "Do not place constraints on dN/dS classes.", 190 "Split", "Assign all dN/dS classes to negative/neutral/positive clases."); 191 if (stratDNDS<0) 192 { 193 return 0; 194 } 195 } 196 else 197 { 198 stratDNDS = 0; 199 } 200 201 if (stratDNDS) 202 { 203 respM = -1; 204 respN = 0; 205 respP = 0; 206 207 while (respM<0 || respM > resp) 208 { 209 fprintf (stdout,"Number of negative rate classes [0-",resp,"]:"); 210 fscanf (stdin, "Number", respM); 211 } 212 213 if (respM < resp) 214 { 215 respN = -1; 216 while (respN<0 || respN > resp-respM) 217 { 218 fprintf (stdout,"Number of neutral rate classes [0-",resp-respM,"]:"); 219 fscanf (stdin, "Number", respN); 220 } 221 } 222 223 respP = resp-respN-respM; 224 225 fprintf (stdout, "\nUsing\n\t", respM, " negatively selected classes\n\t", respN, " neutrally evolving classes\n\t", 226 respP, " positively selected classes\n\n"); 227 } 228 229 categDef1 = ""; 230 categDef1 * 1024; 231 232 lfDef = {}; 233 234 for (fileID = 0; fileID < fileCount; fileID += 1) 235 { 236 lfDef[fileID] = ""; 237 lfDef[fileID] * 1024; 238 lfDef[fileID] * "Log("; 239 } 240 241 242 global S_0 := 1.0; 243 S_0:>0.0000001; 244 global NS_0 = 0.1; 245 246 for (mi=1; mi<resp; mi += 1) 247 { 248 categDef1 * ("global S_"+mi+"=0.5;S_"+mi+":>0.0000001;\nglobal NS_"+mi+";\n"); 249 if (randomizeInitValues) 250 { 251 categDef1*("global P_"+mi+" = Random(0.05,0.95);\nP_"+mi+":<1;\n"); 252 categDef1*("global S_"+mi+" = Random(0.05,1);"); 253 } 254 else 255 { 256 categDef1*("global P_"+mi+" = 1/"+(resp+1-mi)+";\nP_"+mi+":<1;\n"); 257 } 258 } 259 260 freqStrMx = {resp,1}; 261 if (resp>1) 262 { 263 freqStrMx[0] = "P_1"; 264 265 for (mi=1; mi<resp-1; mi=mi+1) 266 { 267 freqStrMx[mi] = ""; 268 for (mi2=1;mi2<=mi;mi2=mi2+1) 269 { 270 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi2+")"; 271 } 272 freqStrMx[mi] = freqStrMx[mi]+"P_"+(mi+1); 273 } 274 freqStrMx[mi] = ""; 275 276 for (mi2=1;mi2<mi;mi2=mi2+1) 277 { 278 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi2+")"; 279 } 280 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi+")"; 281 } 282 else 283 { 284 freqStrMx[0] = "1"; 285 } 286 287 categDef1*( "\n\nglobal c_scale:=S_0*" + freqStrMx[0]); 288 289 for (mi=1; mi<resp; mi=mi+1) 290 { 291 categDef1*( "+S_"+mi+"*" + freqStrMx[mi]); 292 } 293 294 for (mi=0; mi<resp; mi=mi+1) 295 { 296 for (fileID = 0; fileID < fileCount; fileID = fileID + 1) 297 { 298 if (mi) 299 { 300 lfDef[fileID] * "+"; 301 } 302 lfDef[fileID]*(freqStrMx[mi]+"*SITE_LIKELIHOOD["+(fileID*resp+mi)+"]"); 303 } 304 305 } 306 307 categDef1 * ";"; 308 309 for (mi=0; mi<resp; mi+=1) 310 { 311 categDef1 * ("\nglobal R_"+mi+"=1;NS_"+mi+":=R_"+mi+"*S_"+mi+";\n"); 312 } 313 314 if (randomizeInitValues == 2) 315 { 316 enteredValues = {resp,3}; 317 318 fprintf (stdout, "AC="); 319 fscanf (stdin,"Number",thisAlpha); 320 categDef1 * ("AC:="+thisAlpha+";"); 321 fprintf (stdout, "AT="); 322 fscanf (stdin,"Number",thisAlpha); 323 categDef1 * ("AT:="+thisAlpha+";"); 324 fprintf (stdout, "CG="); 325 fscanf (stdin,"Number",thisAlpha); 326 categDef1 * ("CG:="+thisAlpha+";"); 327 fprintf (stdout, "CT="); 328 fscanf (stdin,"Number",thisAlpha); 329 categDef1 * ("CT:="+thisAlpha+";"); 330 fprintf (stdout, "GT="); 331 fscanf (stdin,"Number",thisAlpha); 332 categDef1 * ("GT:="+thisAlpha+";"); 333 334 for (mi = 0; mi < resp; mi = mi+1) 335 { 336 fprintf (stdout, "Alpha for rate class ", mi+1, ":"); 337 fscanf (stdin,"Number",thisAlpha); 338 fprintf (stdout, "Beta for rate class ", mi+1, ":"); 339 fscanf (stdin,"Number",thisBeta); 340 fprintf (stdout, "Weight for rate class ", mi+1, ":"); 341 fscanf (stdin,"Number",thisP); 342 enteredValues [mi][0] = thisAlpha; 343 enteredValues [mi][1] = thisBeta; 344 enteredValues [mi][2] = thisP; 345 } 346 currentMult = 1-enteredValues[0][2]; 347 for (mi = 1; mi < resp-1; mi = mi+1) 348 { 349 enteredValues[mi][2] = enteredValues[mi][2]/currentMult; 350 currentMult = currentMult*(1-enteredValues[mi][2]); 351 } 352 353 for (mi = 0; mi < resp-1; mi = mi+1) 354 { 355 categDef1 * ("S_"+mi+":>0;S_" + mi + ":=" + enteredValues[mi][0] + ";NS_" +mi+ ":=" + enteredValues[mi][1] + ";P_" + (mi+1) + ":=" +enteredValues[mi][2]+";"); 356 } 357 categDef1 * ("S_"+mi+":>0;S_"+mi+":="+enteredValues[mi][0]+";NS_"+mi+":="+enteredValues[mi][1]+";"); 358 } 359 360 if (stratDNDS) 361 { 362 for (mi=respM; mi<respN+respM; mi=mi+1) 363 { 364 categDef1 * ("\nR_"+mi+":=1;"); 365 } 366 367 if (randomizeInitValues) 368 { 369 for (mi=0; mi<respM; mi=mi+1) 370 { 371 categDef1 * ("\nR_"+mi+":<1;R_"+mi+"="+Random(0.05,0.95)+";"); 372 } 373 374 375 for (mi=respM+respN; mi<resp; mi=mi+1) 376 { 377 categDef1 * ("\nR_"+mi+":>1;R_"+mi+"="+Random(1.05,10)+";"); 378 } 379 } 380 else 381 { 382 for (mi=0; mi<respM; mi=mi+1) 383 { 384 categDef1 * ("\nR_"+mi+":<1;R_"+mi+"=1/"+(1+mi)+";"); 385 } 386 387 for (mi=respM+respN; mi<resp; mi=mi+1) 388 { 389 categDef1 * ("\nR_"+mi+":>1;R_"+mi+"=2+"+(mi-respM-respN)+";"); 390 } 391 } 392 } 393 else 394 { 395 if (randomizeInitValues == 1) 396 { 397 for (mi=0; mi<resp; mi=mi+1) 398 { 399 categDef1 * ("\nR_"+mi+"="+Random(0.05,1.75)+";"); 400 } 401 } 402 else 403 { 404 if (randomizeInitValues == 0) 405 { 406 for (mi=0; mi<resp; mi=mi+1) 407 { 408 categDef1 * ("\nR_"+mi+"="+(0.1+0.3*mi)+";"); 409 } 410 } 411 } 412 } 413 414 categDef1 * 0; 415 lfDef1 = ""; 416 lfDef1 * 128; 417 lfDef1 * "\""; 418 419 for (fileID = 0; fileID < fileCount; fileID = fileID + 1) 420 { 421 lfDef[fileID] * ")"; 422 lfDef[fileID] * 0; 423 lfDef1 * lfDef[fileID]; 424 if (fileID < fileCount - 1) 425 { 426 lfDef1 * "+"; 427 } 428 } 429 430 lfDef1 * "\""; 431 lfDef1 * 0; 432 433 ExecuteCommands (categDef1); 434 435 ModelMatrixDimension = 64; 436 for (h = 0 ;h<64; h=h+1) 437 { 438 if (_Genetic_Code[h]==10) 439 { 440 ModelMatrixDimension = ModelMatrixDimension-1; 441 } 442 } 443 444 445 if (randomizeInitValues < 2) 446 { 447 SetDialogPrompt ("Save resulting fit to:"); 448 fprintf (PROMPT_FOR_FILE, CLEAR_FILE); 449 resToPath = LAST_FILE_PATH; 450 } 451 452 if (doCFFreqs == 1) 453 { 454 455 function BuildCodonFrequencies (obsF) 456 { 457 PIStop = 1.0; 458 result = {ModelMatrixDimension,1}; 459 hshift = 0; 460 461 for (h=0; h<64; h=h+1) 462 { 463 first = h$16; 464 second = h%16$4; 465 third = h%4; 466 if (_Genetic_Code[h]==10) 467 { 468 hshift = hshift+1; 469 PIStop = PIStop-obsF[first][0]*obsF[second][1]*obsF[third][2]; 470 continue; 471 } 472 result[h-hshift][0]=obsF[first][0]*obsF[second][1]*obsF[third][2]; 473 } 474 return result*(1.0/PIStop); 475 } 476 paramFreqs = observedFreq; 477 vectorOfFrequencies = BuildCodonFrequencies(paramFreqs); 478 479 } 480 else 481 { 482 ExecuteAFile ("TemplateModels/MGFreqsEstimator.ibf"); 483 BuildCodonFrequencies (paramFreqs, "vectorOfFrequencies"); 484 } 485 486 ExecuteCommands (nucModelString+"\nModel nucModel = (nucModelMatrix,observedFreqSingle);"); 487 488 if (randomizeInitValues == 2) 489 { 490 fileID = ReportDistributionString(resp); 491 fprintf (stdout, "\nUsing the following rate distribution\n", fileID, "\n"); 492 } 493 494 nlfDef = ""; 495 nlfDef * 128; 496 nlfDef * "LikelihoodFunction nuc_lf = ("; 497 498 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 499 { 500 ExecuteCommands ("DataSetFilter nucFilter_"+fileID+" = CreateFilter (filteredData_"+fileID+",1);"); 501 ExecuteCommands ("Tree nucTree_"+fileID+" = treeStringNoModels_"+fileID+";"); 502 if (fileID > 1) 503 { 504 nlfDef * ","; 505 } 506 nlfDef * ("nucFilter_"+fileID+",nucTree_"+fileID); 507 } 508 nlfDef * 0; 509 510 ExecuteCommands (nlfDef + ");"); 511 Optimize (nuc_res, nuc_lf); 512 513 GetString (nucExpStr, nucModel, -1); 514 515 516 lfParts = ""; 517 lfParts * 128; 518 lfParts * "LikelihoodFunction lf = (filteredData_1,tree_1_0"; 519 520 521 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 522 { 523 for (part = 0; part < resp; part = part + 1) 524 { 525 if (fileID == 1) 526 { 527 ExecuteCommands ("PopulateModelMatrix(\"rate_matrix_"+part+"\",paramFreqs,\"S_"+part+"/c_scale\",\"NS_"+part+"/c_scale\",aaRateMultipliers);"); 528 ExecuteCommands ("Model MG94MODEL_"+part+"= (rate_matrix_"+part+",vectorOfFrequencies,0);"); 529 if (part == 0) 530 { 531 GetString (codonExpStr, MG94MODEL_0, -1); 532 synRate = 1; 533 t = 1; 534 //saveS = S_0; 535 //S_0 = 1; 536 saveNS = NS_0; 537 NS_0 = 0.25; 538 GetString (c_scale_constr,c_scale,-2); 539 c_scale = 1; 540 nucL = Eval (nucExpStr); 541 codonL = Eval (codonExpStr); 542 NS_0 = saveNS; 543 //S_0 = saveS; 544 ExecuteCommands ("c_scale:=" + c_scale_constr + ";"); 545 global codonFactor = nucL/codonL; 546 547 } 548 for (_modelID = 1; _modelID < bivariateFitHasMultipleCladeRates; _modelID += 1) 549 { 550 _cladeRate = "clade_" + _modelID + "_NS_" + part; 551 ExecuteCommands ("global " + _cladeRate + " = 1;"); 552 ExecuteCommands ("PopulateModelMatrix(\"rate_matrix_0_clade" + _modelID + "_" + part+"\",paramFreqs,\"S_"+part+"/c_scale\",\"" + _cladeRate + "*NS_"+part+"/c_scale\",aaRateMultipliers);"); 553 ExecuteCommands ("Model MG94MODEL_" + part+"_CLADE_" + _modelID + " = (rate_matrix_0_clade" + _modelID + "_" + part+",vectorOfFrequencies,0);"); 554 } 555 } 556 else 557 { 558 ExecuteCommands ("UseModel (MG94MODEL_"+part+");"); 559 } 560 561 treeID = "tree_"+fileID+"_"+part; 562 ExecuteCommands ("Tree "+treeID+"=treeString_" + fileID +";"); 563 if (branchLengths) 564 { 565 ExecuteCommands ("ReplicateConstraint (\"this1.?.synRate:=this2.?.t__/codonFactor\","+treeID+",nucTree_"+fileID+");"); 566 } 567 else 568 { 569 if (part == 0) 570 { 571 ExecuteCommands ("bnames = BranchName(nucTree_" + fileID + ",-1);"); 572 nlfDef * 128; 573 for (lc = 0; lc < Columns (bnames); lc=lc+1) 574 { 575 nlfDef * (treeID + "." + bnames[lc] + ".synRate = nucTree_" + fileID + "." + bnames[lc] + ".t/codonFactor;"); 576 } 577 nlfDef * 0; 578 ExecuteCommands (nlfDef); 579 } 580 else 581 { 582 ExecuteCommands ("ReplicateConstraint (\"this1.?.synRate:=this2.?.synRate\","+treeID+",tree_1_0);"); 583 } 584 } 585 if (part || fileID > 1) 586 { 587 lfParts = lfParts + ",filteredData_" + fileID + "," + treeID; 588 } 589 } 590 } 591 592 lfParts * 0; 593 ExecuteCommands (lfParts + "," + lfDef1 + ");"); 594 595 USE_LAST_RESULTS = 1; 596 597} 598else 599{ 600 SetDialogPrompt ("Choose a model fit:"); 601 ExecuteAFile (PROMPT_FOR_FILE); 602 603 LF_NEXUS_EXPORT_EXTRA = "bivariateFitHasMultipleCladeRates = " + bivariateFitHasMultipleCladeRates + ";"; 604 605 GetInformation(vars,"^P_[0-9]+$"); 606 resp = Columns (vars)+1; 607 USE_LAST_RESULTS = 1; 608 SKIP_CONJUGATE_GRADIENT = 1; 609 610 ChoiceList (runKind, "Optimization type:", 1, SKIP_NONE, 611 "Unconstrained", "No constraints on rate classes", 612 "Constant dS", "Synonymous rates are constant", 613 "No positive selection", "All dS <= dN", 614 "HKY85", "Nucleotide substitution rates are forced to conform to HKY85"); 615 616 GetString (lfInfo, lf, -1); 617 fileCount = Columns (lfInfo["Datafilters"])/resp; 618 619 GetInformation (branchLengths, "^codonFactor$"); 620 621 //Sfprintf (stdout, "\n\n",branchLengths,"\n\n"); 622 623 branchLengths = (Columns (branchLengths)>0); 624 625 /* do we have codonFactor? */ 626 627 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 628 { 629 ExecuteCommands ("totalCodonCount=totalCodonCount+filteredData_"+fileID+".sites;totalCharCount=totalCharCount+filteredData_"+fileID+".sites*filteredData_"+fileID+".species;"); 630 } 631 632 633 if (runKind < 0) 634 { 635 return 0; 636 } 637 638 639 if (runKind == 0) 640 { 641 resToPath = LAST_FILE_PATH; 642 } 643 else 644 { 645 if (runKind == 1) 646 { 647 resToPath = LAST_FILE_PATH + ".dN_only"; 648 649 c_scale := 1; 650 for (k = 0; k <resp; k=k+1) 651 { 652 ExecuteCommands("global R_"+k+"= Max(NS_"+k+",0.00001)/Max(S_"+k+",0.00001);"); 653 ExecuteCommands("S_"+k+":=1;"); 654 ExecuteCommands("NS_"+k+"=R_"+k+"*S_"+k+";"); 655 } 656 } 657 else 658 { 659 if (runKind == 2) 660 { 661 resToPath = LAST_FILE_PATH + ".NN"; 662 663 for (k = 0; k <resp; k=k+1) 664 { 665 ExecuteCommands("global R_"+k+";R_"+k+":<1;R_"+k+"=Max(NS_"+k+",0.00001)/Max(S_"+k+",0.00001);"); 666 ExecuteCommands("NS_"+k+":=R_"+k+"*S_"+k+";"); 667 } 668 } 669 else 670 { 671 resToPath = LAST_FILE_PATH + ".HKY85"; 672 673 mi = (1+CT)/2; 674 mi2 = ; 675 676 AC = (AC+AT+CG+GT)/(2*(1+CT)); 677 CT := 1; 678 AT := AC; 679 CG := AC; 680 GT := AC; 681 } 682 } 683 } 684 685 freqStrMx = {resp,1}; 686 freqStrMx[0] = "P_1"; 687 688 for (mi=1; mi<resp-1; mi=mi+1) 689 { 690 freqStrMx[mi] = ""; 691 for (mi2=1;mi2<=mi;mi2=mi2+1) 692 { 693 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi2+")"; 694 } 695 freqStrMx[mi] = freqStrMx[mi]+"P_"+(mi+1); 696 } 697 freqStrMx[mi] = ""; 698 699 for (mi2=1;mi2<mi;mi2=mi2+1) 700 { 701 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi2+")"; 702 } 703 freqStrMx[mi] = freqStrMx[mi]+"(1-P_"+mi+")"; 704} 705 706Optimize (res,lf); 707SKIP_CONJUGATE_GRADIENT= 0; 708USE_LAST_RESULTS = 0; 709 710LIKELIHOOD_FUNCTION_OUTPUT = 2; 711sumPath = resToPath + ".summary"; 712 713LogL = res[1][0]; 714 715/* check for branch length approximator */ 716 717GetString (paramList, lf, -1); 718 719degF = Columns(paramList["Global Independent"]) 720 - branchLengths; /* remove one more for codon scaling factor, if nuc branch lengths are used */ 721 722if (doCFFreqs == 1) 723{ 724 degF += 9; 725} 726 727for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 728{ 729 degF += Eval ("BranchCount(tree_" + fileID + "_0)") + Eval ("TipCount(tree_" + fileID + "_0)"); 730} 731 732 733AIC = 2*(degF-LogL); 734AICc = 2*(degF*totalCharCount/(totalCharCount-degF-1) - LogL); 735 736fprintf (stdout, "\n\nModel fit summary\n", 737 "\nLog likelihood:", Format (LogL, 15, 5), 738 "\nParameters :", Format (degF, 15, 0), 739 "\nAIC :", Format (AIC, 15, 5), 740 "\nc-AIC :", Format (AICc, 15, 5),"\n" 741); 742 743fprintf (sumPath, CLEAR_FILE,"Model fit summary\n", 744 "\nLog likelihood:", Format (LogL, 15, 5), 745 "\nParameters :", Format (degF, 15, 0), 746 "\nAIC :", Format (AIC, 15, 5), 747 "\nc-AIC :", Format (AICc, 15, 5),"\n"); 748 749fileID = ReportDistributionString(resp); 750 751bivariateReturnAVL = {}; 752bivariateReturnAVL ["LogL"] = LogL; 753bivariateReturnAVL ["Parameters"] = degF; 754bivariateReturnAVL ["AIC"] = AIC; 755bivariateReturnAVL ["cAIC"] = AICc; 756bivariateReturnAVL ["Rates"] = resp; 757bivariateReturnAVL ["AC"] = AC; 758bivariateReturnAVL ["AT"] = AT; 759bivariateReturnAVL ["CG"] = CG; 760bivariateReturnAVL ["CT"] = CT; 761bivariateReturnAVL ["GT"] = GT; 762bivariateReturnAVL ["Estimates"] = reportMx; 763 764 765if (randomizeInitValues != 2) 766{ 767 fprintf (stdout, fileID); 768 fprintf (sumPath, fileID); 769 770 LIKELIHOOD_FUNCTION_OUTPUT = 7; 771 if (runType == 1) 772 { 773 fprintf (resToPath,CLEAR_FILE,lf); 774 } 775 else 776 { 777 fprintf (resToPath,lf); 778 } 779 780 /*sumPath = resToPath + ".posterior"; 781 ConstructCategoryMatrix (catMat, lf, COMPLETE); 782 783 posteriorProbs = {totalCodonCount, resp}; 784 785 catMatS = ""; 786 catMatS * 1024; 787 788 for (mi2 = 0; mi2 < resp; mi2=mi2+1) 789 { 790 if (mi2) 791 { 792 catMatS * ","; 793 } 794 catMatS * ("Rate class " + (mi2+1) + " omega = " + reportMx[mi2][2]); 795 } 796 797 siteOffset = 0; 798 for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 799 { 800 ExecuteCommands ("thisFilterSize=filteredData_" + fileID + ".sites;"); 801 for (mi=0; mi<thisFilterSize; mi=mi+1) 802 { 803 rs = 0; 804 for (mi2 = 0; mi2 < resp; mi2=mi2+1) 805 { 806 posteriorProbs[mi+siteOffset][mi2] = catMat[mi+siteOffset+mi2*totalCodonCount] * reportMx[mi2][3]; 807 rs = rs + posteriorProbs[mi+siteOffset][mi2]; 808 } 809 catMatS * "\n"; 810 for (mi2 = 0; mi2 < resp; mi2=mi2+1) 811 { 812 posteriorProbs[mi+siteOffset][mi2] = posteriorProbs[mi+siteOffset][mi2]/rs; 813 if (mi2) 814 { 815 catMatS * ","; 816 } 817 catMatS * (""+posteriorProbs[mi+siteOffset][mi2]); 818 } 819 } 820 siteOffset = siteOffset + thisFilterSize; 821 } 822 catMatS * 0; 823 fprintf (sumPath, CLEAR_FILE, catMatS);*/ 824} 825 826 827return bivariateReturnAVL; 828