binaries/src/clustalw/src/tree/ClusterTree.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 #ifdef HAVE_CONFIG_H
   7     #include "config.h"
   8 #endif
   9 #include "ClusterTree.h"
  10 #include "../general/utils.h"
  11 #include "dayhoff.h"
  12 #include "RandomGenerator.h"
  13 #include <math.h>
  14 #include <sstream>
  15 #include "../general/OutputFile.h"
  16
  17 namespace clustalw
  18 {
  19
  20 ClusterTree::ClusterTree()
  21  : numSeqs(0),
  22    firstSeq(0),
  23    lastSeq(0)
  24 {
  25     bootstrapPrompt = "\nEnter name for bootstrap output file  ";
  26     bootstrapFileTypeMsg = "Bootstrap output";
  27 }
  28
  29
  30
  31 /** ****************************************************************************************
  32  *                          The rest are private functions.                                *
  33  *                                                                                         *
  34  *                                                                                         *
  35  *                                                                                         *
  36  *******************************************************************************************/
  37
  38
  39
  40 void ClusterTree::distanceMatrixOutput(ofstream* outFile, clustalw::DistMatrix* matToPrint,
  41                                         clustalw::Alignment *alignPtr)
  42 {
  43     if(outFile == NULL || !outFile->is_open())
  44     {
  45         clustalw::utilityObject->error("Cannot output the distance matrix, file is not open\n");
  46         return;
  47     }
  48
  49     int i, j;
  50     int _maxNames = alignPtr->getMaxNames();
  51     (*outFile) << setw(6) << lastSeq - firstSeq + 1;
  52
  53     for (i = 1; i <= lastSeq - firstSeq + 1; i++)
  54     {
  55         (*outFile) << "\n" << left << setw(_maxNames) << alignPtr->getName(i) << " ";
  56         for (j = 1; j <= lastSeq - firstSeq + 1; j++)
  57         {
  58             (*outFile) << " " << setw(6) << setprecision(3) << fixed << (*matToPrint)(i, j);
  59             if (j % 8 == 0)
  60             {
  61                 if (j != lastSeq - firstSeq + 1)
  62                 {
  63                     (*outFile) << "\n";
  64                 }
  65                 if (j != lastSeq - firstSeq + 1)
  66                 {
  67                     (*outFile) << "          ";
  68                 }
  69             }
  70         }
  71     }
  72 }
  73
  74 void ClusterTree::overspillMessage(int overspill,int totalDists)
  75 {
  76     std::stringstream ssOverSpill;
  77     ssOverSpill << overspill << " of the distances out of a total of " << totalDists
  78          << "\n were out of range for the distance correction." << "\n"
  79          << "\n SUGGESTIONS: 1) remove the most distant sequences"
  80          << "\n           or 2) use the PHYLIP package"
  81          << "\n           or 3) turn off the correction."
  82          << "\n Note: Use option 3 with caution! With this degree"
  83          << "\n of divergence you will have great difficulty"
  84          << "\n getting robust and reliable trees." << "\n\n";
  85     string message;
  86     ssOverSpill >> message;
  87     clustalw::utilityObject->warning(message.c_str());
  88 }
  89
  90 /*
  91  * The function treeGapDelete flags all the positions that have a gap in any sequence.
  92  *
  93  */
  94 // NOTE there is something wrong with using _lenFirstSeq. But this is what old clustal does.
  95 void ClusterTree::treeGapDelete(clustalw::Alignment *alignPtr)
  96 {
  97     int seqn;
  98     int posn;
  99     int _maxAlnLength = alignPtr->getMaxAlnLength();
 100     int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
 101     int _gapPos1 = clustalw::userParameters->getGapPos1();
 102     int _gapPos2 = clustalw::userParameters->getGapPos2();
 103
 104     treeGaps.resize(_maxAlnLength + 1);
 105
 106     for (posn = 1; posn <= _lenFirstSeq; ++posn)
 107     {
 108         treeGaps[posn] = 0;
 109         for (seqn = 1; seqn <= lastSeq - firstSeq + 1; ++seqn)
 110         {
 111             const vector<int>* _seqM = alignPtr->getSequence(seqn + firstSeq - 1);
 112
 113             if(posn > alignPtr->getSeqLength(seqn + firstSeq - 1))
 114             {
 115                 break; // Dont read locations that cannot be read!
 116             }
 117             if (((*_seqM)[posn] == _gapPos1) || ((*_seqM)[posn] == _gapPos2))
 118             {
 119                 treeGaps[posn] = 1;
 120                 break;
 121             }
 122         }
 123     }
 124 }
 125
 126 int ClusterTree::dnaDistanceMatrix(ofstream* treeFile, clustalw::Alignment *alignPtr)
 127 {
 128     int m, n;
 129     int j, i;
 130     int res1, res2;
 131     int overspill = 0;
 132     double p, q, e, a, b, k;
 133
 134     treeGapDelete(alignPtr); // flag positions with gaps (tree_gaps[i] = 1 )
 135
 136     if (verbose)
 137     {
 138         (*treeFile) << "\n";
 139         (*treeFile) <<  "\n DIST   = percentage divergence (/100)";
 140         (*treeFile) << "\n p      = rate of transition (A <-> G; C <-> T)";
 141         (*treeFile) << "\n q      = rate of transversion";
 142         (*treeFile) << "\n Length = number of sites used in comparison";
 143         (*treeFile) << "\n";
 144         if (clustalw::userParameters->getTossGaps())
 145         {
 146             (*treeFile) << "\n All sites with gaps (in any sequence) deleted!";
 147             (*treeFile) << "\n";
 148         }
 149         if (clustalw::userParameters->getKimura())
 150         {
 151             (*treeFile) << "\n Distances corrected by Kimura's 2 parameter model:";
 152             (*treeFile) << "\n\n Kimura, M. (1980)";
 153             (*treeFile) << " A simple method for estimating evolutionary ";
 154             (*treeFile) << "rates of base";
 155             (*treeFile) << "\n substitutions through comparative studies of ";
 156             (*treeFile) << "nucleotide sequences.";
 157             (*treeFile) << "\n J. Mol. Evol., 16, 111-120.";
 158             (*treeFile) << "\n\n";
 159         }
 160     }
 161
 162     int _numSeqs = alignPtr->getNumSeqs();
 163     quickDistMat.reset(new clustalw::DistMatrix(_numSeqs + 1));
 164     int _lenFirstSeq = alignPtr->getSeqLength(firstSeq);
 165     int _gapPos1 = clustalw::userParameters->getGapPos1();
 166     int _gapPos2 = clustalw::userParameters->getGapPos2();
 167     int lenSeqM, lenSeqN;
 168      // for every pair of sequence
 169     for (m = 1; m < lastSeq - firstSeq + 1; ++m)
 170     {
 171         const vector<int>* _seqM = alignPtr->getSequence(m + firstSeq - 1);
 172         lenSeqM = alignPtr->getSeqLength(m + firstSeq - 1);
 173
 174         for (n = m + 1; n <= lastSeq - firstSeq + 1; ++n)
 175         {
 176             const vector<int>* _seqN = alignPtr->getSequence(n + firstSeq - 1);
 177             lenSeqN = alignPtr->getSeqLength(n + firstSeq - 1);
 178
 179             p = q = e = 0.0;
 180             quickDistMat->SetAt(m, n, 0.0);
 181             quickDistMat->SetAt(n ,m, 0.0);
 182
 183             for (i = 1; i <= _lenFirstSeq; ++i)
 184             {
 185                 j = bootPositions[i];
 186                 if (clustalw::userParameters->getTossGaps() && (treeGaps[j] > 0))
 187                 {
 188                     goto skip;
 189                 }
 190                 // gap position
 191 /** *******************************************************************************
 192  * BUG!!!!!!!      NOTE this was found for protein. Presuming the same here       *
 193  * NOTE: the following if statements were coded in so as to produce               *
 194  * the same distance results as the old clustal. Old clustal compares             *
 195  * up to the length of the first sequence. If this is longer than the             *
 196  * other sequences, then the -3 and 0's are compared at the end of the            *
 197  * array. These should not be compared, but I need to stick to this to            *
 198  * produce the same results as the old version for testing!                       *
 199  **********************************************************************************/
 200                 if(j > lenSeqM)
 201                 {
 202                     if(j == lenSeqM + 1)
 203                     {
 204                         res1 = -3;
 205                     }
 206                     else
 207                     {
 208                         res1 = 0;
 209                     }
 210                 }
 211                 else
 212                 {
 213                     res1 = (*_seqM)[j];
 214                 }
 215                 if(j > lenSeqN)
 216                 {
 217                     if(j == lenSeqN + 1)
 218                     {
 219                         res2 = -3;
 220                     }
 221                     else
 222                     {
 223                         res2 = 0;
 224                     }
 225                 }
 226                 else
 227                 {
 228                     res2 = (*_seqN)[j];
 229                 }
 230                 if ((res1 == _gapPos1) || (res1 == _gapPos2) || (res2 == _gapPos1)
 231                     || (res2 == _gapPos2))
 232                 {
 233                     goto skip;
 234                 }
 235                 // gap in a seq
 236                 if (!clustalw::userParameters->getUseAmbiguities())
 237                 {
 238                     if (isAmbiguity(res1) || isAmbiguity(res2))
 239                     {
 240                         goto skip;
 241                     }
 242                 }
 243                 // ambiguity code in a seq
 244                 e = e+1.0;
 245                 if (res1 != res2)
 246                 {
 247                     if (transition(res1, res2))
 248                     {
 249                         p = p + 1.0;
 250                     }
 251                     else
 252                     {
 253                         q = q + 1.0;
 254                     }
 255                 }
 256                 skip: ;
 257             }
 258
 259
 260             // Kimura's 2 parameter correction for multiple substitutions
 261
 262             if (!clustalw::userParameters->getKimura())
 263             {
 264                 if (e == 0)
 265                 {
 266                     cerr << "\n WARNING: sequences " << m << " and " << n
 267                          << " are non-overlapping\n";
 268                     k = 0.0;
 269                     p = 0.0;
 270                     q = 0.0;
 271                 }
 272                 else
 273                 {
 274                     k = (p + q) / e;
 275                     if (p > 0.0)
 276                     {
 277                         p = p / e;
 278                     }
 279                     else
 280                     {
 281                         p = 0.0;
 282                     }
 283                     if (q > 0.0)
 284                     {
 285                         q = q / e;
 286                     }
 287                     else
 288                     {
 289                         q = 0.0;
 290                     }
 291                 }
 292                 quickDistMat->SetAt(m, n, k);
 293                 quickDistMat->SetAt(n ,m, k);
 294                 if (verbose)
 295                 {
 296                     (*treeFile) << setw(4) << m << " vs." << setw(4) << n << ":  DIST =  "
 297                                 << setw(4) << fixed << setprecision(4)
 298                                 << k << "; p = " << fixed << setprecision(4) << p << "; q = "
 299                                 << fixed << setprecision(4) << q << "; length = " << setw(6)
 300                                 << fixed << setprecision(0) << e << "\n";
 301                 }
 302             }
 303             else
 304             {
 305                 if (e == 0)
 306                 {
 307                     cerr << "\n WARNING: sequences " << m << " and " << n
 308                          << " are non-overlapping\n";;
 309                     p = 0.0;
 310                     q = 0.0;
 311                 }
 312                 else
 313                 {
 314                     if (p > 0.0)
 315                     {
 316                         p = p / e;
 317                     }
 318                     else
 319                     {
 320                         p = 0.0;
 321                     }
 322                     if (q > 0.0)
 323                     {
 324                         q = q / e;
 325                     }
 326                     else
 327                     {
 328                         q = 0.0;
 329                     }
 330                 }
 331
 332                 if (((2.0 *p) + q) == 1.0)
 333                 {
 334                     a = 0.0;
 335                 }
 336                 else
 337                 {
 338                     a = 1.0 / (1.0 - (2.0 *p) - q);
 339                 }
 340
 341                 if (q == 0.5)
 342                 {
 343                     b = 0.0;
 344                 }
 345                 else
 346                 {
 347                     b = 1.0 / (1.0 - (2.0 *q));
 348                 }
 349
 350                 // watch for values going off the scale for the correction.
 351                 if ((a <= 0.0) || (b <= 0.0))
 352                 {
 353                     overspill++;
 354                     k = 3.5; // arbitrary high score
 355                 }
 356                 else
 357                 {
 358                     k = 0.5 * log(a) + 0.25 * log(b);
 359                 }
 360                 quickDistMat->SetAt(m, n, k);
 361                 quickDistMat->SetAt(n ,m, k);
 362                 if (verbose)
 363                 // if screen output
 364                 {
 365                     (*treeFile) << setw(4) << m << " vs." << setw(4) << n
 366                                 << ":  DIST =  " << fixed << setprecision(4)
 367                                 << k << "; p = " << fixed <<  setprecision(4) << p << "; q = "
 368                                 << fixed << setprecision(4) << q << "; length = " << setw(6)
 369                                 << fixed << setprecision(0) << e << "\n";
 370                 }
 371
 372             }
 373         }
 374     }
 375     return overspill; // return the number of off-scale values
 376 }
 377
 378 int ClusterTree::protDistanceMatrix(ofstream* treeFile, clustalw::Alignment *alignPtr)
 379 {
 380     int m, n;
 381     int j, i;
 382     int res1, res2;
 383     int overspill = 0;
 384     double p, e, k, tableEntry;
 385
 386     treeGapDelete(alignPtr); // flag positions with gaps (tree_gaps[i] = 1 )
 387
 388     if (verbose)
 389     {
 390         (*treeFile) <<  "\n";
 391         (*treeFile) << "\n DIST   = percentage divergence (/100)";
 392         (*treeFile) << "\n Length = number of sites used in comparison";
 393         (*treeFile) << "\n\n";
 394         if (clustalw::userParameters->getTossGaps())
 395         {
 396             (*treeFile) << "\n All sites with gaps (in any sequence) deleted";
 397             (*treeFile) << "\n";
 398         }
 399         if (clustalw::userParameters->getKimura())
 400         {
 401             (*treeFile) <<
 402                 "\n Distances up to 0.75 corrected by Kimura's empirical method:";
 403             (*treeFile) << "\n\n Kimura, M. (1983)";
 404             (*treeFile) << " The Neutral Theory of Molecular Evolution.";
 405             (*treeFile) <<
 406                 "\n Page 75. Cambridge University Press, Cambridge, England.";
 407             (*treeFile) << "\n\n";
 408         }
 409     }
 410     int _numSeqs = alignPtr->getNumSeqs();
 411     int _lenSeq1 = alignPtr->getSeqLength(1);
 412     quickDistMat.reset(new clustalw::DistMatrix(_numSeqs + 1));
 413     int _gapPos1 = clustalw::userParameters->getGapPos1();
 414     int _gapPos2 = clustalw::userParameters->getGapPos2();
 415     int lenSeqM, lenSeqN;
 416     // for every pair of sequence
 417     for (m = 1; m < _numSeqs; ++m)
 418     {
 419         const vector<int>* _seqM = alignPtr->getSequence(m);
 420         lenSeqM = alignPtr->getSeqLength(m);
 421         for (n = m + 1; n <= _numSeqs; ++n)
 422         {
 423             const vector<int>* _seqN = alignPtr->getSequence(n);
 424             lenSeqN = alignPtr->getSeqLength(n);
 425             p = e = 0.0;
 426             quickDistMat->SetAt(m, n, 0.0);
 427             quickDistMat->SetAt(n ,m, 0.0);
 428             for (i = 1; i <= _lenSeq1; ++i) // It may be this here!
 429             {
 430                 j = bootPositions[i];
 431                 if (clustalw::userParameters->getTossGaps() && (treeGaps[j] > 0))
 432                 {
 433                     goto skip;
 434                 }
 435                 // gap position
 436 /** *******************************************************************************
 437  * BUG!!!!!!!                                                                     *
 438  * NOTE: the following if statements were coded in so as to produce               *
 439  * the same distance results as the old clustal. Old clustal compares             *
 440  * up to the length of the first sequence. If this is longer than the             *
 441  * other sequences, then the -3 and 0's are compared at the end of the            *
 442  * array. These should not be compared, but I need to stick to this to            *
 443  * produce the same results as the old version for testing!                       *
 444  **********************************************************************************/
 445                 if(j > lenSeqM)
 446                 {
 447                     if(j == lenSeqM + 1)
 448                     {
 449                         res1 = -3;
 450                     }
 451                     else
 452                     {
 453                         res1 = 0;
 454                     }
 455                 }
 456                 else
 457                 {
 458                     res1 = (*_seqM)[j];
 459                 }
 460                 if(j > lenSeqN)
 461                 {
 462                     if(j == lenSeqN + 1)
 463                     {
 464                         res2 = -3;
 465                     }
 466                     else
 467                     {
 468                         res2 = 0;
 469                     }
 470                 }
 471                 else
 472                 {
 473                     res2 = (*_seqN)[j];
 474                 }
 475                 if ((res1 == _gapPos1) || (res1 == _gapPos2) || (res2 == _gapPos1)
 476                     || (res2 == _gapPos2))
 477                 {
 478                     goto skip;
 479                 }
 480                 // gap in a seq
 481                 e = e + 1.0;
 482                 if (res1 != res2)
 483                 {
 484                     p = p + 1.0;
 485                 }
 486                 skip: ;
 487             }
 488
 489             if (p <= 0.0)
 490             {
 491                 k = 0.0;
 492             }
 493             else
 494             {
 495                 k = p / e;
 496             }
 497
 498             if (clustalw::userParameters->getKimura())
 499             {
 500                 if (k < 0.75)
 501                 {
 502                     // use Kimura's formula
 503                     if (k > 0.0)
 504                     {
 505                         k =  - log(1.0 - k - (k * k / 5.0));
 506                     }
 507                 }
 508                 else
 509                 {
 510                     if (k > 0.930)
 511                     {
 512                         overspill++;
 513                         k = 10.0; // arbitrarily set to 1000%
 514                     }
 515                     else // dayhoff_pams is from dayhoff.h file
 516                     {
 517                         tableEntry = (k * 1000.0) - 750.0;
 518                         k = (double)dayhoff_pams[(int)tableEntry];
 519                         k = k / 100.0;
 520                     }
 521                 }
 522             }
 523
 524             quickDistMat->SetAt(m, n, k);
 525             quickDistMat->SetAt(n ,m, k);
 526             if (verbose)
 527             {
 528                 (*treeFile) << setw(4) << m << " vs." << setw(4) << n
 529                             << "  DIST = " << fixed << setprecision(4)
 530                             << k << ";  length = " << setw(6)
 531                             << setprecision(0) << e << "\n";
 532             }
 533         }
 534     }
 535     return overspill;
 536 }
 537
 538 bool ClusterTree::isAmbiguity(int c)
 539 {
 540     int i;
 541     char codes[] = "ACGTU";
 542
 543     if (clustalw::userParameters->getUseAmbiguities() == true)
 544     {
 545         return false;
 546     }
 547
 548     for (i = 0; i < 5; i++)
 549         if (clustalw::userParameters->getAminoAcidCode(c) == codes[i])
 550         {
 551             return false;
 552         }
 553
 554     return true;
 555 }
 556
 557 /*
 558  * The function calcPercIdentity calculates the percent identity of the sequences
 559  * and outputs it to a the file pfile. NOTE this is not used at the moment. It was in
 560  * the old code, but there was no way to access it from the menu. This may change.
 561  */
 562 void ClusterTree::calcPercIdentity(ofstream* pfile, clustalw::Alignment *alignPtr)
 563 {
 564     clustalw::DistMatrix percentMat;
 565
 566     float ident;
 567     int nmatch;
 568
 569     int val1, val2;
 570
 571     int i,j,k, length_longest;
 572     int length_shortest;
 573
 574     int rs = 0, rl = 0;
 575     // findout sequence length, longest and shortest;
 576     length_longest = 0;
 577     length_shortest = 0;
 578
 579     int _numSeqs = alignPtr->getNumSeqs();
 580     int _seqLength;
 581     for (i = 1; i <= _numSeqs; i++)
 582     {
 583         _seqLength = alignPtr->getSeqLength(i);
 584         if (length_longest < _seqLength)
 585         {
 586             length_longest = _seqLength;
 587             rs = i;
 588         }
 589         if (length_shortest > _seqLength)
 590         {
 591             length_shortest = _seqLength;
 592             rl = i;
 593         }
 594     }
 595
 596     percentMat.ResizeRect(_numSeqs + 1);
 597     nmatch = 0;
 598     int _lenSeqI, _lenSeqJ;
 599     int _maxAA = clustalw::userParameters->getMaxAA();
 600
 601     for (i = 1; i <= _numSeqs; i++)
 602     {
 603         const vector<int>* _seqI = alignPtr->getSequence(i);
 604         _lenSeqI = alignPtr->getSeqLength(i);
 605
 606         for (j = i; j <= _numSeqs; j++)
 607         {
 608             const vector<int>* _seqJ = alignPtr->getSequence(j);
 609             _lenSeqJ = alignPtr->getSeqLength(j);
 610
 611             cout << "\n           " << alignPtr->getName(j) << " ";
 612             ident = 0;
 613             nmatch = 0;
 614             for(k = 1; k <= length_longest; k++)
 615             {
 616                 if((k > _lenSeqI) || (k > _lenSeqJ))
 617                 {
 618                     break;
 619                 }
 620                 val1 = (*_seqI)[k];
 621                 val2 = (*_seqJ)[k];
 622
 623                 if (((val1 < 0) || (val1 > _maxAA)) || ((val2 < 0) || (val2 > _maxAA)))
 624                 {
 625                     continue; // residue = '-';
 626                 }
 627                 if (val1 == val2)
 628                 {
 629                     ident++ ;
 630                     nmatch++;
 631                 }
 632                 else
 633                 {
 634                     nmatch++ ;
 635                 }
 636             }
 637             ident = ident/nmatch * 100.0 ;
 638             percentMat.SetAt(i, j, ident);
 639             percentMat.SetAt(j, i, ident);
 640         }
 641     }
 642
 643     int _maxNameSize = alignPtr->getMaxNames();
 644
 645     (*pfile) << "#\n#\n#  Percent Identity  Matrix - created by Clustal"
 646              << clustalw::userParameters->getRevisionLevel() << " \n#\n#\n";
 647     for(i = 1; i <= _numSeqs; i++)
 648     {
 649         (*pfile) << "\n " << right << setw(5) << i << ": ";
 650         (*pfile) << left << setw(_maxNameSize) << alignPtr->getName(i);
 651
 652         for(j = 1; j <= _numSeqs; j++)
 653         {
 654             (*pfile) << setw(8) << right << fixed << setprecision(0) << percentMat(i, j);
 655         }
 656     }
 657     (*pfile) << "\n";
 658
 659 }
 660
 661 void ClusterTree::compareTree(PhyloTree* tree1, PhyloTree* tree2, vector<int>* hits, int n)
 662 {
 663     int i, j, k;
 664     int nhits1, nhits2;
 665
 666     for (i = 1; i <= n - 3; i++)
 667     {
 668         for (j = 1; j <= n - 3; j++)
 669         {
 670             nhits1 = 0;
 671             nhits2 = 0;
 672             for (k = 1; k <= n; k++)
 673             {
 674                 if (tree1->treeDesc[i][k] == tree2->treeDesc[j][k])
 675                 {
 676                     nhits1++;
 677                 }
 678                 if (tree1->treeDesc[i][k] != tree2->treeDesc[j][k])
 679                 {
 680                     nhits2++;
 681                 }
 682             }
 683             if ((nhits1 == lastSeq - firstSeq + 1) || (nhits2 == lastSeq -
 684                 firstSeq + 1))
 685             {
 686                 (*hits)[i]++;
 687             }
 688         }
 689     }
 690 }
 691
 692 /**
 693  * NOTE this will go into the OutputFile class and will not be needed here anymore.
 694  *
 695  */
 696 /*string ClusterTree::getOutputFileName(const string prompt, string path,
 697                                       const string fileExtension)
 698 {
 699     string temp;
 700     string _fileName; // Will return this name.
 701     string message;
 702     _fileName = path + fileExtension;
 703
 704     if(_fileName.compare(clustalw::userParameters->getSeqName()) == 0)
 705     {
 706         cout << "Output file name is the same as input file.\n";
 707         if (clustalw::userParameters->getMenuFlag())
 708         {
 709             message = "\n\nEnter new name to avoid overwriting  [" + _fileName + "]: ";
 710             clustalw::utilityObject->getStr(message, temp);
 711             if(temp != "")
 712             {
 713                 _fileName = temp;
 714             }
 715         }
 716     }
 717     else if (clustalw::userParameters->getMenuFlag())
 718     {
 719
 720         message = prompt + " [" + _fileName + "]";
 721         clustalw::utilityObject->getStr(message, temp);
 722         if(temp != "")
 723         {
 724             _fileName = temp;
 725         }
 726     }
 727     return _fileName;
 728
 729 }*/
 730
 731 bool ClusterTree::transition(int base1, int base2)
 732 {
 733 // assumes that the bases of DNA sequences have been translated as
 734 // a,A = 0;   c,C = 1;   g,G = 2;   t,T,u,U = 3;  N = 4;
 735 // a,A = 0;   c,C = 2;   g,G = 6;   t,T,u,U =17;
 736 // A <--> G  and  T <--> C  are transitions;  all others are transversions.
 737     if (((base1 == 0) && (base2 == 6)) || ((base1 == 6) && (base2 == 0)))
 738     {
 739         return true;
 740     }
 741     // A <--> G
 742     if (((base1 == 17) && (base2 == 2)) || ((base1 == 2) && (base2 == 17)))
 743     {
 744         return true;
 745     }
 746     // T <--> C
 747     return false;
 748 }
 749
 750 /**
 751  * This function is used to open all the bootstrap tree files. It opens them with the
 752  * correct message prompt.
 753  */
 754 bool ClusterTree::openFilesForBootstrap(clustalw::OutputFile* clustalFile, clustalw::OutputFile* phylipFile,
 755                          clustalw::OutputFile* nexusFile, clustalw::TreeNames* treeNames, string* path)
 756 {
 757     if (clustalw::userParameters->getOutputTreeClustal())
 758     {
 759         if(!clustalFile || !clustalFile->openFile(&(treeNames->clustalName),
 760                                   bootstrapPrompt, path, "njb", bootstrapFileTypeMsg))
 761         {
 762             return false;
 763         }
 764     }
 765
 766     if (clustalw::userParameters->getOutputTreePhylip())
 767     {
 768         if(!phylipFile || !phylipFile->openFile(&(treeNames->phylipName),
 769                                 bootstrapPrompt, path, "phb", bootstrapFileTypeMsg))
 770         {
 771             return false;
 772         }
 773     }
 774
 775     if (clustalw::userParameters->getOutputTreeNexus())
 776     {
 777         if(!nexusFile || !nexusFile->openFile(&(treeNames->nexusName),
 778                                 bootstrapPrompt, path, "treb", bootstrapFileTypeMsg))
 779         {
 780             return false;
 781         }
 782     }
 783     return true;
 784 }
 785
 786 bool ClusterTree::openFilesForTreeFromAlignment(clustalw::OutputFile* clustalFile,
 787     clustalw::OutputFile* phylipFile, clustalw::OutputFile* distFile, clustalw::OutputFile* nexusFile, clustalw::OutputFile* pimFile,
 788                             clustalw::TreeNames* treeNames, string* path)
 789 {
 790     if (clustalw::userParameters->getOutputTreeClustal())
 791     {
 792         if(!clustalFile || !clustalFile->openFile(&(treeNames->clustalName),
 793                                   "\nEnter name for CLUSTAL    tree output file  ",
 794                                   path, "nj", "Phylogenetic tree"))
 795         {
 796             return false;
 797         }
 798     }
 799
 800     if (clustalw::userParameters->getOutputTreePhylip())
 801     {
 802         if(!phylipFile || !phylipFile->openFile(&(treeNames->phylipName),
 803                              "\nEnter name for PHYLIP     tree output file  ", path, "ph",
 804                              "Phylogenetic tree"))
 805         {
 806             return false;
 807         }
 808     }
 809
 810     if (clustalw::userParameters->getOutputTreeDistances())
 811     {
 812         if(!distFile || !distFile->openFile(&(treeNames->distName),
 813                                "\nEnter name for distance matrix output file  ",
 814                                path, "dst", "Distance matrix"))
 815         {
 816             return false;
 817         }
 818     }
 819
 820     if (clustalw::userParameters->getOutputTreeNexus())
 821     {
 822         if(!nexusFile || !nexusFile->openFile(&(treeNames->nexusName),
 823                                 "\nEnter name for NEXUS tree output file  ", path,
 824                                      "tre", "Nexus tree"))
 825         {
 826             return false;
 827         }
 828     }
 829
 830     if (clustalw::userParameters->getOutputPim())
 831     {
 832         if(!pimFile || !pimFile->openFile(&(treeNames->pimName),
 833                             "\nEnter name for % Identity matrix output file  ", path, "pim",
 834                             "perc identity"))
 835         {
 836             return false;
 837         }
 838     }
 839     return true;
 840 }
 841
 842 int ClusterTree::calcQuickDistMatForAll(ofstream* clustalFile, ofstream* phylipFile,
 843                                   ofstream* nexusFile, ofstream* pimFile, ofstream* distFile,
 844                                   clustalw::Alignment* alignPtr)
 845 {
 846     int overspill = 0;
 847     bool _DNAFlag = clustalw::userParameters->getDNAFlag();
 848
 849     overspill = calcQuickDistMatForSubSet(clustalFile, phylipFile, nexusFile, alignPtr);
 850
 851     if (pimFile && clustalw::userParameters->getOutputPim())
 852     {
 853         verbose = false; // Turn off file output
 854         if (_DNAFlag)
 855         {
 856             calcPercIdentity(pimFile, alignPtr);
 857         }
 858         else
 859         {
 860             calcPercIdentity(pimFile, alignPtr);
 861         }
 862     }
 863
 864
 865     if (distFile && clustalw::userParameters->getOutputTreeDistances())
 866     {
 867         verbose = false; // Turn off file output
 868         if (_DNAFlag)
 869         {
 870             overspill = dnaDistanceMatrix(distFile, alignPtr);
 871         }
 872         else
 873         {
 874             overspill = protDistanceMatrix(distFile, alignPtr);
 875         }
 876         distanceMatrixOutput(distFile, quickDistMat.get(),
 877                              alignPtr);
 878     }
 879     return overspill;
 880 }
 881
 882 int ClusterTree::calcQuickDistMatForSubSet(ofstream* clustalFile, ofstream* phylipFile,
 883                                           ofstream* nexusFile, clustalw::Alignment* alignPtr,
 884                                           bool inBootLoop)
 885 {
 886     int overspill = 0;
 887     bool _DNAFlag = clustalw::userParameters->getDNAFlag();
 888
 889     if (clustalFile && clustalw::userParameters->getOutputTreeClustal())
 890     {
 891         if(!inBootLoop)
 892         {
 893             verbose = true; // Turn on file output
 894         }
 895         else
 896         {
 897             verbose = false; // Turn off when we are in the loop in bootstrap!
 898         }
 899         if (_DNAFlag)
 900         {
 901             overspill = dnaDistanceMatrix(clustalFile, alignPtr);
 902         }
 903         else
 904         {
 905             overspill = protDistanceMatrix(clustalFile, alignPtr);
 906         }
 907     }
 908
 909     if (phylipFile && clustalw::userParameters->getOutputTreePhylip())
 910     {
 911         verbose = false; // Turn off file output
 912         if (_DNAFlag)
 913         {
 914             overspill = dnaDistanceMatrix(phylipFile, alignPtr);
 915         }
 916         else
 917         {
 918             overspill = protDistanceMatrix(phylipFile, alignPtr);
 919         }
 920     }
 921
 922     if (nexusFile && clustalw::userParameters->getOutputTreeNexus())
 923     {
 924         verbose = false; // Turn off file output
 925         if (_DNAFlag)
 926         {
 927             overspill = dnaDistanceMatrix(nexusFile, alignPtr);
 928         }
 929         else
 930         {
 931             overspill = protDistanceMatrix(nexusFile, alignPtr);
 932         }
 933     }
 934     return overspill;
 935 }
 936
 937 void ClusterTree::printBootstrapHeaderToClustalFile(ofstream* clustalFile)
 938 {
 939     if(clustalFile)
 940     {
 941         (*clustalFile) << "\n\n\t\t\tBootstrap Confidence Limits\n\n";
 942         (*clustalFile) << "\n Random number generator seed = "
 943                        << setw(7)
 944                        << clustalw::userParameters->getBootRanSeed() << "\n";
 945         (*clustalFile) << "\n Number of bootstrap trials   = " << setw(7)
 946                        << clustalw::userParameters->getBootNumTrials() << "\n";
 947         (*clustalFile) << "\n\n Diagrammatic representation of the above tree: \n";
 948         (*clustalFile) << "\n Each row represents 1 tree cycle;";
 949         (*clustalFile) << " defining 2 groups.\n";
 950         (*clustalFile) << "\n Each column is 1 sequence; ";
 951         (*clustalFile) << "the stars in each line show 1 group; ";
 952         (*clustalFile) << "\n the dots show the other\n";
 953         (*clustalFile) << "\n Numbers show occurences in bootstrap samples.";
 954     }
 955 }
 956
 957 void ClusterTree::promptForBoolSeedAndNumTrials()
 958 {
 959     if (clustalw::userParameters->getMenuFlag())
 960     {
 961         unsigned int tempSeed;
 962         tempSeed = clustalw::utilityObject->getInt(
 963                 "\n\nEnter seed no. for random number generator ", 1, 1000,
 964         clustalw::userParameters->getBootRanSeed());
 965         clustalw::userParameters->setBootRanSeed(tempSeed);
 966
 967         clustalw::userParameters->setBootNumTrials(
 968                      clustalw::utilityObject->getInt("\n\nEnter number of bootstrap trials ",
 969                         1, 10000, clustalw::userParameters->getBootNumTrials()));
 970     }
 971 }
 972
 973 void ClusterTree::printErrorMessageForBootstrap(int totalOverspill, int totalDists,
 974                                                 int nfails)
 975 {
 976     cerr << "\n";
 977     cerr << "\n WARNING: " << totalOverspill
 978          << " of the distances out of a total of "
 979          << totalDists << " times" << clustalw::userParameters->getBootNumTrials();
 980     cerr << "\n were out of range for the distance correction.";
 981     cerr << "\n This affected " << nfails << " out of "
 982          << clustalw::userParameters->getBootNumTrials() << " bootstrap trials.";
 983     cerr << "\n This may not be fatal but you have been warned!" << "\n";
 984     cerr << "\n SUGGESTIONS: 1) turn off the correction";
 985     cerr << "\n           or 2) remove the most distant sequences";
 986     cerr << "\n           or 3) use the PHYLIP package.\n\n";
 987
 988     if (clustalw::userParameters->getMenuFlag())
 989     {
 990         string dummy;
 991         clustalw::utilityObject->getStr(string("Press [RETURN] to continue"), dummy);
 992     }
 993 }
 994
 995 bool ClusterTree::checkIfConditionsMet(int numSeqs, int min)
 996 {
 997     if (clustalw::userParameters->getEmpty())
 998     {
 999         clustalw::utilityObject->error("You must load an alignment first");
1000         return false;
1001     }
1002
1003     if (numSeqs < min)
1004     {
1005         clustalw::utilityObject->error("Alignment has only %d sequences", numSeqs);
1006         return false;
1007     }
1008
1009     return true;
1010 }
1011
1012 }