1#include "binomial.ibf"; 2 3/*___________________________________________________________________________________________________________*/ 4 5function PadString (padLength,padChar) 6{ 7 padString = ""; 8 d=padString*padLength; 9 for (padCounter=0;padCounter<padLength;padCounter=padCounter+1) 10 { 11 d=padString*padChar; 12 } 13 d=padString*0; 14 fprintf (stdout,padString); 15 return padString; 16} 17 18/*___________________________________________________________________________________________________________*/ 19 20function PrintASCIITable (dataMatrix, titleMatrix) 21{ 22 columnWidths = {1,Columns(titleMatrix)}; 23 for (counter1=0; counter1<Columns(titleMatrix); counter1 = counter1+1) 24 { 25 counter2 = Abs (titleMatrix[0][counter1])+2; 26 if (counter2<12) 27 { 28 counter2 = 12; 29 } 30 columnWidths[0][counter1] = counter2; 31 } 32 fprintf (stdout, "\n"); 33 for (counter2 = 0; counter2 < Columns (titleMatrix); counter2 = counter2+1) 34 { 35 fprintf (stdout,"+-"); 36 dummy = PadString (columnWidths[0][counter2],"-"); 37 fprintf (stdout,"-"); 38 } 39 fprintf (stdout,"+\n| "); 40 41 for (counter1=0; counter1<Columns(titleMatrix); counter1 = counter1+1) 42 { 43 fprintf (stdout, titleMatrix[counter1]); 44 dummy = PadString (columnWidths[0][counter1]-Abs(titleMatrix[counter1])," "); 45 fprintf (stdout, " | "); 46 } 47 48 fprintf (stdout, "\n"); 49 50 for (counter1=-1; counter1<Rows(dataMatrix); counter1 = counter1 + 1) 51 { 52 if (counter1>=0) 53 { 54 fprintf (stdout,"| "); 55 fprintf (stdout,Format(counter1+1,columnWidths[0][0],0)); 56 for (counter2 = 1; counter2 < Columns (titleMatrix); counter2 = counter2+1) 57 { 58 fprintf (stdout," | "); 59 fprintf (stdout,Format(dataMatrix[counter1][counter2-1],columnWidths[0][counter2],-1)); 60 } 61 fprintf (stdout," "); 62 fprintf (stdout, "|\n"); 63 } 64 for (counter2 = 0; counter2 < Columns (titleMatrix); counter2 = counter2+1) 65 { 66 fprintf (stdout,"+-"); 67 dummy = PadString (columnWidths[0][counter2],"-"); 68 fprintf (stdout,"-"); 69 } 70 fprintf (stdout, "+\n"); 71 } 72 73 return 1; 74} 75 76/*___________________________________________________________________________________________________________*/ 77 78nucCharacters="ACGT"; 79function codonString (ccode) 80{ 81 return nucCharacters[ccode$16]+nucCharacters[(ccode%16)$4]+nucCharacters[ccode%4]; 82} 83 84/*___________________________________________________________________________________________________________*/ 85 86function PrintTableToFile (dataMatrix, titleMatrix, promptOrNot) 87{ 88 SetDialogPrompt ("Export tab separated data to:"); 89 90 if (promptOrNot) 91 { 92 fprintf (PROMPT_FOR_FILE, CLEAR_FILE); 93 } 94 95 bufferString = ""; 96 bufferString * 8192; 97 bufferString * titleMatrix[0][0]; 98 99 for (counter1=1; counter1<Columns(titleMatrix); counter1 = counter1+1) 100 { 101 bufferString * ("\t"+titleMatrix[0][counter1]); 102 } 103 104 for (counter1=0; counter1<Rows(dataMatrix); counter1 = counter1 + 1) 105 { 106 bufferString * ("\n"+dataMatrix[counter1][0]); 107 for (counter2 = 1; counter2 < Columns (dataMatrix); counter2 = counter2+1) 108 { 109 bufferString * ("\t"+dataMatrix[counter1][counter2]); 110 } 111 } 112 113 bufferString * 0; 114 fprintf (LAST_FILE_PATH,bufferString); 115 bufferString = 0; 116 117 return 1; 118} 119 120/*----------------------------------------------------------------------------*/ 121 122ChoiceList (ambChoice, "Treatment of Ambiguities",1,SKIP_NONE, 123 "Averaged","All possible resolutions are considered and averaged.", 124 "Resolved","The most frequent (for that site) resolution is chosen."); 125 126if (ambChoice<0) 127{ 128 return; 129} 130 131if (useCustomCountingBias) 132{ 133 ExecuteAFile("Distances/CodonToolsMain.def"); 134} 135else 136{ 137 ExecuteAFile("Distances/CodonTools.def"); 138} 139 140observedCEFV = {64,1}; 141for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 142{ 143 ExecuteCommands ("HarvestFrequencies (tp, filteredData_"+fileID+",3,3,0);cfs = filteredData_"+fileID+".sites;"); 144 observedCEFV = observedCEFV + tp*(cfs/totalCodonCount); 145} 146 147 148seqToBranchMap = {stateCharCount,1}; 149hShift = 0; 150for (k=0; k<64; k=k+1) 151{ 152 if (_Genetic_Code[k]==10) 153 { 154 hShift = hShift+1; 155 } 156 else 157 { 158 seqToBranchMap[k-hShift] = observedCEFV[k]; 159 } 160} 161observedCEFV = seqToBranchMap; 162 163vOffset = 0; 164resultMatrix = {totalCodonCount,12}; 165 166treeLengthArray = {}; 167 168for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 169{ 170 ExecuteCommands ("LikelihoodFunction tempLF = (filteredData_"+fileID+",codonTree_"+fileID+");"); 171 DataSet dsA = ReconstructAncestors (tempLF); 172 ExecuteCommands ("DataSet dsJoint = Combine(dsA,ds_"+fileID+");"); 173 ExecuteCommands ("DataSetFilter filteredData = CreateFilter (ds_"+fileID+",3,\"\",\"\",GeneticCodeExclusions);"); 174 175 DataSetFilter filteredDataA = CreateFilter (dsA,3,"","",GeneticCodeExclusions); 176 DataSetFilter filteredDataJ = CreateFilter (dsJoint,3,"","",GeneticCodeExclusions); 177 178 ExecuteCommands ("branchNames = BranchName (codonTree_"+fileID+",-1);"); 179 h = Columns (branchNames); 180 181 seqToBranchMap = {h, 2}; 182 /* maps sequence names to branch order in column 1 183 and the other way around in column 2 */ 184 185 for (k=0; k<filteredData.species; k=k+1) 186 { 187 GetString (seqName, filteredData, k); 188 seqToBranchMap[k][0] = -1; 189 for (v=0; v<h; v=v+1) 190 { 191 if (branchNames[v] % seqName) 192 { 193 seqToBranchMap[k][0] = v; 194 seqToBranchMap[v][1] = k; 195 break; 196 } 197 } 198 } 199 200 seqToBranchMap[filteredData.species][0] = h-1; 201 seqToBranchMap[h-1][1] = filteredData.species; 202 203 for (k=1; k<filteredDataA.species; k=k+1) 204 { 205 GetString (seqName, filteredDataA, k); 206 seqToBranchMap[filteredData.species+k][0] = -1; 207 for (v=0; v<h; v=v+1) 208 { 209 if (branchNames[v] % seqName) 210 { 211 seqToBranchMap[k+filteredData.species][0] = v; 212 seqToBranchMap[v][1] = k+filteredData.species; 213 break; 214 } 215 } 216 } 217 218 219 /* total tree length */ 220 221 totalTreeLength = 0; 222 223 ExecuteCommands ("branchLengths = BranchLength(nucTree_"+fileID+",-1);"); 224 for (v=Columns(branchLengths)-1; v>=0; v=v-1) 225 { 226 totalTreeLength = totalTreeLength + branchLengths[v]; 227 } 228 229 treeLengthArray [fileID] = totalTreeLength; 230 231 /* get codon matrix */ 232 233 filteredData.unique_sites = Rows(filteredData.site_freqs) * Columns (fitleredData.site_freqs); 234 235 codonInfo = {filteredData.species, filteredData.unique_sites}; 236 codonInfo2 = {filteredDataA.species, filteredDataA.unique_sites}; 237 238 GetDataInfo (dupInfo, filteredData); 239 GetDataInfo (dupInfoA, filteredDataA); 240 241 matrixTrick = {1,stateCharCount}; 242 matrixTrick2 = {1,stateCharCount}; 243 matrixTrick = matrixTrick["_MATRIX_ELEMENT_COLUMN_"]; 244 matrixTrick2 = matrixTrick2["1"]; 245 246 for (v=0; v<filteredData.unique_sites;v=v+1) 247 { 248 for (h=0; h<filteredData.species;h=h+1) 249 { 250 GetDataInfo (siteInfo, filteredData, h, v); 251 _SITE_ES_COUNT = matrixTrick2 * siteInfo; 252 if (_SITE_ES_COUNT[0] == 1) 253 { 254 siteInfo = matrixTrick * siteInfo; 255 codonInfo[h][v] = siteInfo[0]; 256 } 257 else 258 { 259 codonInfo[h][v] = -1; 260 } 261 } 262 } 263 264 for (v=0; v<filteredDataA.unique_sites;v=v+1) 265 { 266 for (h=0; h<filteredDataA.species;h=h+1) 267 { 268 GetDataInfo (siteInfo, filteredDataA, h, v); 269 if ((matrixTrick2*siteInfo)[0] > 1) 270 { 271 codonInfo2[h][v] = -1; 272 } 273 else 274 { 275 codonInfo2[h][v] = (matrixTrick * siteInfo)[0]; 276 } 277 } 278 } 279 _SITE_RESULTS = {4,filteredData.sites}; 280 ExecuteCommands ("flatTreeRep = Abs (nucTree_"+fileID+");"); 281 282 perBranchMatrix = {Columns(branchNames),2}; 283 284 for (v=0; v<filteredData.sites;v=v+1) 285 { 286 _SITE_ES_COUNT = {stateCharCount,stateCharCount}; 287 _SITE_EN_COUNT = {stateCharCount,stateCharCount}; 288 _SITE_OS_COUNT = {stateCharCount,stateCharCount}; 289 _SITE_ON_COUNT = {stateCharCount,stateCharCount}; 290 291 /* do internal nodes first */ 292 293 k = filteredData.species+1; 294 295 /* first sequence is always the root */ 296 c1 = dupInfoA[v]; 297 if (codonInfo2[1][c1] >= stateCharCount) 298 { 299 continue; 300 } 301 for (h=1; h<filteredDataA.species; h=h+1) 302 { 303 p1 = seqToBranchMap[k][0]; 304 p2 = flatTreeRep[p1]; 305 p2 = seqToBranchMap[p2][1]-filteredData.species; 306 307 branchFactor = branchLengths[p1]/totalTreeLength; 308 309 cd1 = codonInfo2[h] [c1]; 310 cd2 = codonInfo2[p2][c1]; 311 312 if (cd1 >= 0 && cd2 >= 0) 313 { 314 _SITE_OS_COUNT[cd1][cd2] = _SITE_OS_COUNT[cd1][cd2] + 1; 315 _SITE_ON_COUNT[cd1][cd2] = _SITE_ON_COUNT[cd1][cd2] + 1; 316 _SITE_ES_COUNT[cd1][cd2] = _SITE_ES_COUNT[cd1][cd2] + branchFactor; 317 _SITE_EN_COUNT[cd1][cd2] = _SITE_EN_COUNT[cd1][cd2] + branchFactor; 318 319 if (_PAIRWISE_S_[cd1][cd2]) 320 { 321 perBranchMatrix[p1][0] = perBranchMatrix[p1][0]+_OBSERVED_S_[cd1][cd2]/_PAIRWISE_S_[cd1][cd2]; 322 } 323 perBranchMatrix[p1][1] = perBranchMatrix[p1][1]+_OBSERVED_NS_[cd1][cd2]/_PAIRWISE_NS_[cd1][cd2]; 324 } 325 326 k=k+1; 327 } 328 329 /* now do the leaves */ 330 331 observedCEFV = {{0}}; 332 333 c2 = dupInfo[v]; 334 for (h=0; h<filteredData.species; h=h+1) 335 { 336 p1 = seqToBranchMap[h][0]; 337 p2 = flatTreeRep[p1]; 338 p2 = seqToBranchMap[p2][1]-filteredData.species; 339 340 cd2 = codonInfo2[p2][c1]; 341 cd1 = codonInfo[h] [c2]; 342 343 branchFactor = branchLengths[p1]/totalTreeLength; 344 345 if (cd1>=0) 346 /* no ambiguities */ 347 { 348 /*p2 = flatTreeRep[p1]; 349 fprintf (stdout, "Change from ", codonString (cd1), " to ", codonString (cd2), " along ", branchNames[p1], " and ", branchNames[p2], "(",branchFactor,")\n"); 350 */ 351 _SITE_OS_COUNT[cd1][cd2] = _SITE_OS_COUNT[cd1][cd2] + 1; 352 _SITE_ON_COUNT[cd1][cd2] = _SITE_ON_COUNT[cd1][cd2] + 1; 353 _SITE_ES_COUNT[cd1][cd2] = _SITE_ES_COUNT[cd1][cd2] + branchFactor; 354 _SITE_EN_COUNT[cd1][cd2] = _SITE_EN_COUNT[cd1][cd2] + branchFactor; 355 356 if (_PAIRWISE_S_[cd1][cd2]) 357 { 358 perBranchMatrix[p1][0] = perBranchMatrix[p1][0]+_OBSERVED_S_[cd1][cd2]/_PAIRWISE_S_[cd1][cd2]; 359 } 360 perBranchMatrix[p1][1] = perBranchMatrix[p1][1]+_OBSERVED_NS_[cd1][cd2]/_PAIRWISE_NS_[cd1][cd2]; 361 } 362 else 363 /* ambiguities here */ 364 { 365 GetDataInfo (ambInfo, filteredData, h, c2); 366 if (Rows(observedCEFV) == 1) 367 { 368 siteFilter = ""+(v*3)+"-"+(v*3+2); 369 DataSetFilter filteredDataSite = CreateFilter (dsJoint,3,siteFilter,"",GeneticCodeExclusions); 370 HarvestFrequencies (observedCEFV,filteredDataSite,3,3,0); 371 tempMx = {stateCharCount,1}; 372 373 hShift = 0; 374 375 for (k=0; k<64; k=k+1) 376 { 377 if (_Genetic_Code[k]==10) 378 { 379 hShift = hShift+1; 380 } 381 else 382 { 383 tempMx[k-hShift] = observedCEFV[k]; 384 } 385 } 386 observedCEFV = tempMx; 387 } 388 389 weightFactor = matrixTrick2*ambInfo; 390 if (weightFactor[0]<stateCharCount) 391 { 392 ambInfo = ambInfo$observedCEFV; 393 394 if (ambChoice) 395 { 396 weightFactor = 0; 397 tempMx = -1; 398 for (k=0; k<stateCharCount; k=k+1) 399 { 400 if (ambInfo[k]>weightFactor) 401 { 402 weightFactor = ambInfo[k]; 403 tempMx = k; 404 } 405 } 406 if (tempMx>=0) 407 { 408 _SITE_OS_COUNT[tempMx][cd2] = _SITE_OS_COUNT[tempMx][cd2] + 1; 409 _SITE_ON_COUNT[tempMx][cd2] = _SITE_ON_COUNT[tempMx][cd2] + 1 ; 410 _SITE_ES_COUNT[tempMx][cd2] = _SITE_ES_COUNT[tempMx][cd2] + branchFactor; 411 _SITE_EN_COUNT[tempMx][cd2] = _SITE_EN_COUNT[tempMx][cd2] + branchFactor; 412 if (_PAIRWISE_S_[tempMx][cd2]) 413 { 414 perBranchMatrix[p1][0] = perBranchMatrix[p1][0]+_OBSERVED_S_[tempMx][cd2]/_PAIRWISE_S_[tempMx][cd2]; 415 } 416 perBranchMatrix[p1][1] = perBranchMatrix[p1][1]+_OBSERVED_NS_[tempMx][cd2]/_PAIRWISE_NS_[tempMx][cd2]; 417 } 418 } 419 else 420 { 421 weightFactor = matrixTrick2*ambInfo; 422 weightFactor = weightFactor[0]; 423 if (weightFactor) 424 { 425 ambInfo = ambInfo * (1/weightFactor); 426 for (k=0; k<stateCharCount; k=k+1) 427 { 428 weightFactor = ambInfo[k]; 429 if (weightFactor>0) 430 { 431 _SITE_OS_COUNT[k][cd2] = _SITE_OS_COUNT[k][cd2] + weightFactor; 432 _SITE_ON_COUNT[k][cd2] = _SITE_ON_COUNT[k][cd2] + weightFactor; 433 _SITE_ES_COUNT[k][cd2] = _SITE_ES_COUNT[k][cd2] + weightFactor*branchFactor; 434 _SITE_EN_COUNT[k][cd2] = _SITE_EN_COUNT[k][cd2] + weightFactor*branchFactor; 435 436 if (_PAIRWISE_S_[k][cd2]) 437 { 438 perBranchMatrix[p1][0] = perBranchMatrix[p1][0]+weightFactor*_OBSERVED_S_[k][cd2]/_PAIRWISE_S_[k][cd2]; 439 } 440 perBranchMatrix[p1][1] = perBranchMatrix[p1][1]+weightFactor*_OBSERVED_NS_[k][cd2]/_PAIRWISE_NS_[k][cd2]; 441 } 442 } 443 } 444 } 445 } 446 } 447 } 448 449 _SITE_OS_COUNT = matrixTrick2*(_OBSERVED_S_$_SITE_OS_COUNT)*Transpose(matrixTrick2); 450 _SITE_ON_COUNT = matrixTrick2*(_OBSERVED_NS_$_SITE_ON_COUNT)*Transpose(matrixTrick2); 451 _SITE_ES_COUNT = matrixTrick2*(_PAIRWISE_S_$_SITE_ES_COUNT)*Transpose(matrixTrick2); 452 _SITE_EN_COUNT = matrixTrick2*(_PAIRWISE_NS_$_SITE_EN_COUNT)*Transpose(matrixTrick2); 453 454 shiftedV = v+vOffset; 455 456 resultMatrix[shiftedV][0] = _SITE_OS_COUNT[0]; 457 resultMatrix[shiftedV][1] = _SITE_ON_COUNT[0]; 458 resultMatrix[shiftedV][2] = _SITE_ES_COUNT[0]; 459 resultMatrix[shiftedV][3] = _SITE_EN_COUNT[0]; 460 461 weight = _SITE_EN_COUNT[0]+_SITE_ES_COUNT[0]; 462 463 p = _SITE_ES_COUNT[0]/weight; 464 465 resultMatrix[shiftedV][5] = p; 466 467 p2 = resultMatrix[shiftedV][0]+resultMatrix[shiftedV][1]; 468 469 resultMatrix[shiftedV][7] = resultMatrix[shiftedV][1]/resultMatrix[shiftedV][3]; 470 471 if (resultMatrix[shiftedV][2]) 472 { 473 resultMatrix[shiftedV][6] = resultMatrix[shiftedV][0]/resultMatrix[shiftedV][2]; 474 resultMatrix[shiftedV][8] = resultMatrix[shiftedV][7]-resultMatrix[shiftedV][6]; 475 resultMatrix[shiftedV][11] = resultMatrix[shiftedV][8]/totalTreeLength; 476 } 477 478 if (p2) 479 { 480 resultMatrix[shiftedV][4] = resultMatrix[shiftedV][0]/p2; 481 if (distribChoice) 482 { 483 p3 = p2$1; 484 p4 = ((resultMatrix[shiftedV][0]*6)+0.5)$1; 485 486 waP = 0; 487 waN = 0; 488 489 if (p3 != p2) 490 { 491 /* use averaging */ 492 fprintf (stdout, "Interpolating at site ", v, ": ", p3, " ", p3+1, " ", p2, "\n"); 493 494 w1 = p2-p3; 495 w2 = 1-w1; 496 497 for (kk2 = 0; kk2 < stateCharCount; kk2 = kk2+1) 498 { 499 ambInfo = SUPPORT_MATRIX_LIST[kk2][v]; 500 if (ambInfo>0.000001) 501 { 502 subMatrixAVL = simulatedNullAVL[kk2]; 503 subMatrixAVL1 = subMatrixAVL[p3]; 504 subMatrixAVL2 = subMatrixAVL[p3+1]; 505 if (Abs(subMatrixAVL1) && Abs(subMatrixAVL2)) 506 { 507 pv1 = subMatrixAVL1[p4+1][1]; 508 pv2 = subMatrixAVL2[p4+1][1]; 509 if (subMatrixAVL1[0][0] >= 100 && subMatrixAVL2[0][0] >= 100) 510 { 511 waP = waP + (pv1*w1+pv2*w2)*ambInfo; 512 if (p4 == 0) 513 { 514 waN = waN + ambInfo; 515 } 516 else 517 { 518 waN = waN + (1 + w1*(pv1 - 2*subMatrixAVL1[p4][1]) + w2*(pv2 - 2*subMatrixAVL2[p4][1])) * ambInfo; 519 } 520 } 521 else 522 { 523 fprintf (stdout, "Bailing at site ", v, " codon ", kk2, " sim count ", subMatrixAVL[0][0], "\n"); 524 break; 525 } 526 } 527 else 528 { 529 fprintf (stdout, "Bailing at site ", v, " codon ", kk2 , " sim count 0 \n"); 530 break; 531 } 532 } 533 } 534 } 535 else 536 { 537 for (kk2 = 0; kk2 < stateCharCount; kk2 = kk2+1) 538 { 539 ambInfo = SUPPORT_MATRIX_LIST[kk2][v]; 540 if (ambInfo>0.000001) 541 { 542 subMatrixAVL = simulatedNullAVL[kk2]; 543 subMatrixAVL = subMatrixAVL[p3]; 544 if (Abs(subMatrixAVL)) 545 { 546 if (subMatrixAVL[0][0] >= 100) 547 { 548 waP = waP + subMatrixAVL[p4+1][1]*ambInfo; 549 if (p4 == 0) 550 { 551 waN = waN + ambInfo; 552 } 553 else 554 { 555 waN = waN + (1 + subMatrixAVL[p4+1][1] - 2*subMatrixAVL[p4][1]) * ambInfo; 556 } 557 } 558 else 559 { 560 fprintf (stdout, "Bailing at site ", v, " codon ", kk2, " sim count ", subMatrixAVL[0][0], "\n"); 561 break; 562 } 563 } 564 else 565 { 566 fprintf (stdout, "Bailing at site ", v, " codon ", kk2 , " sim count 0 \n"); 567 break; 568 } 569 } 570 } 571 } 572 573 574 575 if (kk2 == stateCharCount) 576 { 577 resultMatrix[shiftedV][9] = waP; 578 resultMatrix[shiftedV][10] = waN; 579 continue; 580 } 581 } 582 resultMatrix[shiftedV][9] = extendedBinTail (p2,p,resultMatrix[shiftedV][0]); 583 if (resultMatrix[shiftedV][0]>=1) 584 { 585 resultMatrix[shiftedV][10] = 1-extendedBinTail(p2,p,resultMatrix[shiftedV][0]-1); 586 } 587 else 588 { 589 resultMatrix[shiftedV][10] = 1-extendedBinTail (p2,p,0); 590 } 591 } 592 } 593 vOffset = vOffset + filteredData.sites; 594} 595 596perBranchMatrix = perBranchMatrix * (1/totalCodonCount); 597labelMatrix = {1,12}; 598labelMatrix[0] = "Observed S Changes"; 599labelMatrix[1] = "Observed NS Changes"; 600labelMatrix[2] = "E[S Sites]"; 601labelMatrix[3] = "E[NS Sites]"; 602labelMatrix[4] = "Observed S. Prop."; 603labelMatrix[5] = "P{S}"; 604labelMatrix[6] = "dS"; 605labelMatrix[7] = "dN"; 606labelMatrix[8] = "dN-dS"; 607labelMatrix[9] = "P{S leq. observed}"; 608labelMatrix[10] = "P{S geq. observed}"; 609labelMatrix[11] = "Scaled dN-dS"; 610 611sigLevel = -1; 612 613for (fileID = 1; fileID <= fileCount; fileID = fileID + 1) 614{ 615 fprintf (stdout, "\nTotal length (subs/nuc/unit time) for tree ", fileID, ": ", Format (treeLengthArray[fileID],10,5)); 616} 617 618fprintf (stdout, "\n"); 619 620while ((sigLevel<=0)||(sigLevel>=1)) 621{ 622 fprintf (stdout, "\nSignificance level for a site to be classified as positively/negatively selected?"); 623 fscanf (stdin, "Number", sigLevel); 624} 625 626posSelected = 0; 627negSelected = 0; 628 629p = Rows(resultMatrix); 630 631for (p2=0; p2<p; p2=p2+1) 632{ 633 v = resultMatrix [p2][8]; 634 if (v>0) 635 { 636 if (resultMatrix [p2][9] < sigLevel) 637 { 638 posSelected = posSelected+1; 639 } 640 } 641 else 642 { 643 if (v<0) 644 { 645 if (resultMatrix [p2][10] < sigLevel) 646 { 647 negSelected = negSelected+1; 648 } 649 } 650 } 651} 652 653selLabelMatrix = {{"Index","Site Index","dN-dS","p-value"}}; 654 655if (posSelected) 656{ 657 psMatrix = {posSelected, 3}; 658 h = 0; 659 for (p2=0; p2<p; p2=p2+1) 660 { 661 v = resultMatrix [p2][8]; 662 if (v>0) 663 { 664 if (resultMatrix [p2][9] < sigLevel) 665 { 666 psMatrix[h][0] = p2+1; 667 psMatrix[h][1] = v; 668 psMatrix[h][2] = resultMatrix [p2][9]; 669 h = h+1; 670 } 671 } 672 } 673 674 fprintf (stdout,"\n******* FOUND ", posSelected, " POSITIVELY SELECTED SITES ********\n\n"); 675 dummy = PrintASCIITable (psMatrix, selLabelMatrix); 676} 677else 678{ 679 fprintf (stdout,"\n******* FOUND NO POSITIVELY SELECTED SITES ********\n\n"); 680} 681 682if (negSelected) 683{ 684 psMatrix = {negSelected, 3}; 685 h = 0; 686 for (p2=0; p2<p; p2=p2+1) 687 { 688 v = resultMatrix [p2][8]; 689 if (v<0) 690 { 691 if (resultMatrix [p2][10] < sigLevel) 692 { 693 psMatrix[h][0] = p2+1; 694 psMatrix[h][1] = v; 695 psMatrix[h][2] = resultMatrix [p2][10]; 696 h = h+1; 697 } 698 } 699 } 700 701 fprintf (stdout,"\n******* FOUND ", negSelected, " NEGATIVELY SELECTED SITES ********\n\n"); 702 dummy = PrintASCIITable (psMatrix, selLabelMatrix); 703} 704else 705{ 706 fprintf (stdout,"\n******* FOUND NO NEGATIVELY SELECTED SITES ********\n\n"); 707} 708 709ChoiceList (outputChoice, "Output Options",1,SKIP_NONE, 710 "ASCII Table", "Output is printed to the console as an ASCII table.", 711 "Export to File","Output is spooled to a tab separated file.", 712 "Chart","A HYPHY chart window is displayed (Mac/Doze Only)."); 713if (outputChoice<0) 714{ 715 return; 716} 717 718if (outputChoice==0) 719{ 720 dummy = PrintASCIITable (resultMatrix, labelMatrix); 721} 722else 723{ 724 if (outputChoice == 1) 725 { 726 SetDialogPrompt ("Save result table to"); 727 dummy = PrintTableToFile (resultMatrix, labelMatrix, !pipeThroughFlag); 728 } 729 else 730 { 731 OpenWindow (CHARTWINDOW,{{"Rates by sites"} 732 {"labelMatrix"}, 733 {"resultMatrix"}, 734 {"Contrast Bars"}, 735 {"Index"}, 736 {labelMatrix[6]+";"+labelMatrix[7]}, 737 {"Site Index"}, 738 {"dN"}, 739 {"dS"}, 740 {"0"}}, 741 "SCREEN_WIDTH-60;SCREEN_HEIGHT-50;30;50"); 742 } 743} 744