binaries/src/clustalw/src/tree/Tree.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 <sstream>
  10 #include "Tree.h"
  11
  12 namespace clustalw
  13 {
  14
  15 /**
  16  *
  17  * @param firstSeq
  18  * @param lastSeq
  19  * @param sweight
  20  */
  21
  22 void Tree::calcSeqWeights(int firstSeq, int lastSeq, vector<int>* sweight)
  23 {
  24     if((int)sweight->size() < lastSeq - 1)
  25     {
  26         sweight->resize(lastSeq - 1);
  27     }
  28
  29     int i, _nSeqs;
  30     int temp, sum;
  31     int *weight;
  32     /*
  33      * If there are more than three sequences....
  34      */
  35     _nSeqs = lastSeq - firstSeq;
  36     if ((_nSeqs >= 2) && (clustalw::userParameters->getDistanceTree() == true) &&
  37         (clustalw::userParameters->getNoWeights() == false))
  38     {
  39         /*
  40          * Calculate sequence weights based on Phylip tree.
  41          */
  42
  43         weight = new int[lastSeq + 1];
  44
  45         for (i = firstSeq; i < lastSeq; i++)
  46         {
  47             weight[i] = calcWeight(i);
  48         }
  49
  50         /*
  51          * Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
  52          */
  53
  54         sum = 0;
  55         for (i = firstSeq; i < lastSeq; i++)
  56         {
  57             sum += weight[i];
  58         }
  59
  60         if (sum == 0)
  61         {
  62             for (i = firstSeq; i < lastSeq; i++)
  63             {
  64                 weight[i] = 1;
  65             }
  66             sum = i;
  67         }
  68
  69         for (i = firstSeq; i < lastSeq; i++)
  70         {
  71             (*sweight)[i] = (weight[i] * INT_SCALE_FACTOR) / sum;
  72             if ((*sweight)[i] < 1)
  73             {
  74                 (*sweight)[i] = 1;
  75             }
  76         }
  77
  78         delete []weight;
  79         weight = NULL;
  80     }
  81
  82     else
  83     {
  84         /*
  85          * Otherwise, use identity weights.
  86          */
  87         temp = INT_SCALE_FACTOR / _nSeqs;
  88         // AW 2009-05-28: goes wrong if we have more than
  89         // INT_SCALE_FACTOR seqs. if so, set to 1, just as above
  90         if (temp < 1)
  91             temp = 1;
  92
  93         for (i = firstSeq; i < lastSeq; i++)
  94         {
  95             (*sweight)[i] = temp;
  96         }
  97     }
  98
  99 }
 100
 101 /**
 102  *
 103  * @param alignPtr
 104  * @param treeFileName
 105  * @param firstSeq
 106  * @param lastSeq
 107  * @return
 108  */
 109 int Tree::readTree(clustalw::Alignment* alignPtr, const string& treeFileName, int firstSeq, int lastSeq)
 110 {
 111
 112     if(alignPtr == NULL || firstSeq < 0 || lastSeq < 1)
 113     {
 114         return 0;
 115     }
 116     char c;
 117     string name1 = "", name2 = "";
 118     int i, j, k;
 119     bool found;
 120
 121     numSeq = 0;
 122     nnodes = 0;
 123     ntotal = 0;
 124     rootedTree = true;
 125
 126     // Need to check what happens if I try to open a file that doesnt exist!
 127     if(treeFileName.empty())
 128     {
 129         return 0; // Cannot open, empty string!
 130     }
 131     else
 132     {
 133         file.open(treeFileName.c_str(), ifstream::in);
 134         if(!file.is_open())
 135         {
 136             clustalw::utilityObject->error("Cannot open output file [%s]\n", treeFileName.c_str());
 137         }
 138     }
 139
 140     skipSpace(&file);
 141     charFromFile = file.get();
 142
 143     if (charFromFile != '(')
 144     {
 145         clustalw::utilityObject->error("Wrong format in tree file %s\n", treeFileName.c_str());
 146         return 0;
 147     }
 148     file.seekg(0, ios::beg); // Put get pointer back to the begining!
 149
 150     clustalw::userParameters->setDistanceTree(true);
 151
 152
 153     // Allocate memory for tree
 154
 155     // Allocate memory.
 156     nptr = new TreeNode*[3 * (lastSeq - firstSeq + 1)];
 157     ptrs = new TreeNode*[3 * (lastSeq - firstSeq + 1)];
 158     lptr = new TreeNode*[lastSeq - firstSeq + 1];
 159     olptr = new TreeNode*[lastSeq + 1];
 160
 161     seqTree = avail();
 162     setInfo(seqTree, NULL, 0, string(""), 0.0);
 163
 164     createTree(seqTree, NULL, &file);
 165     file.close();
 166
 167     if (numSeq != lastSeq - firstSeq)
 168     {
 169         stringstream ss;
 170         ss << "tree not compatible with alignment\n(" << lastSeq - firstSeq
 171         << " sequences in alignment and "<< numSeq <<" in tree\n";
 172         string errorMes;
 173         ss >> errorMes;
 174         clustalw::utilityObject->error(errorMes.c_str());
 175         return 0;
 176     }
 177
 178     // If the tree is unrooted, reroot the tree - ie. minimise the difference
 179     // between the mean root->leaf distances for the left and right branches of
 180     // the tree.
 181
 182     if (clustalw::userParameters->getDistanceTree() == false)
 183     {
 184         if (rootedTree == false)
 185         {
 186             clustalw::utilityObject->error("input tree is unrooted and has no distances.\nCannot align sequences");
 187             return 0;
 188         }
 189     }
 190
 191     if (rootedTree == false)
 192     {
 193         root = reRoot(seqTree, lastSeq - firstSeq + 1);
 194     }
 195     else
 196     {
 197         root = seqTree;
 198     }
 199
 200     // calculate the 'order' of each node.
 201     int nameLength;
 202     string nameSeq;
 203     orderNodes();
 204
 205     if (numSeq >= 2)
 206     {
 207         // If there are more than three sequences....
 208         // assign the sequence nodes (in the same order as in the alignment file)
 209         for (i = firstSeq; i < lastSeq; i++)
 210         {
 211             nameLength = alignPtr->getName(i + 1).length();
 212             nameSeq = alignPtr->getName(i + 1);
 213             name1 = "";
 214             name2 = "";
 215             if (nameLength > clustalw::MAXNAMES)
 216             {
 217                 stringstream ss;
 218                 ss << "name " << nameSeq << " is too long for PHYLIP tree format (max "
 219                    << clustalw::MAXNAMES << " chars)";
 220                 string msg;
 221                 ss >> msg;
 222                 clustalw::utilityObject->warning(msg.c_str());// Mark change 17-5-07
 223             }
 224
 225             for (k = 0; k < nameLength && k < clustalw::MAXNAMES; k++)
 226             {
 227                 c = nameSeq[k];
 228                 if ((c > 0x40) && (c < 0x5b))
 229                 {
 230                     c = c | 0x20;
 231                 }
 232                 if (c == ' ')
 233                 {
 234                     c = '_';
 235                 }
 236                 name2 += c;
 237             }
 238             found = false;
 239
 240             for (j = 0; j < numSeq; j++)
 241             {
 242                 name1 = "";
 243                 for (k = 0; k < (int)lptr[j]->name.length() && k < clustalw::MAXNAMES; k++)
 244                 {
 245                     c = lptr[j]->name[k];
 246                     if ((c > 0x40) && (c < 0x5b))
 247                     {
 248                         c = c | 0x20;
 249                     }
 250
 251                     name1 += c;
 252                 }
 253
 254                 if (name1.compare(name2) == 0)
 255                 {
 256                     olptr[i] = lptr[j];
 257                     found = true;
 258                 }
 259             }
 260
 261             if (found == false)
 262             {
 263                 utilityObject->error("tree not compatible with alignment:\n %s not found\n", name2.c_str());
 264                 return 0;
 265             }
 266         }
 267
 268     }
 269     return (1);
 270 }
 271
 272 /**
 273  *
 274  * @param firstSeq
 275  * @param lastSeq
 276  * @return
 277  */
 278 auto_ptr<AlignmentSteps> Tree::createSets(int firstSeq, int lastSeq)
 279 {
 280     auto_ptr<AlignmentSteps> progAlignSteps;
 281     progAlignSteps.reset(new AlignmentSteps);
 282
 283     int i, j, _nSeqs;
 284
 285     numSets = 0;
 286     _nSeqs = lastSeq - firstSeq;
 287     if (_nSeqs >= 2)
 288     {
 289         // If there are more than three sequences....
 290         groups = new int[_nSeqs + 1];
 291         groupSeqs(root, groups, _nSeqs, progAlignSteps.get());
 292
 293         delete []groups;
 294         groups  = NULL;
 295
 296     }
 297
 298     else
 299     {
 300         groups = new int[_nSeqs + 1];
 301         for (i = 0; i < _nSeqs - 1; i++)
 302         {
 303             for (j = 0; j < _nSeqs; j++)
 304                 if (j <= i)
 305                 {
 306                     groups[j] = 1;
 307                 }
 308                 else if (j == i + 1)
 309                 {
 310                     groups[j] = 2;
 311                 }
 312                 else
 313                 {
 314                     groups[j] = 0;
 315                 }
 316
 317             progAlignSteps->saveSet(_nSeqs, groups);
 318         }
 319         delete []groups;
 320         groups  = NULL;
 321     }
 322
 323     return progAlignSteps;
 324 }
 325
 326 /**
 327  * calcSimilarities changes the distMat.
 328  * @param alignPtr
 329  * @param distMat
 330  * @return
 331  */
 332 int Tree::calcSimilarities(clustalw::Alignment* alignPtr, clustalw::DistMatrix* distMat)
 333 {
 334     int depth = 0, i, j, k, n;
 335     bool found;
 336     int nerrs, seq1[MAXERRS], seq2[MAXERRS];
 337     TreeNode *p,  **pathToRoot;
 338     float dist;
 339     float *distToNode, badDist[MAXERRS];
 340     double **dmat;
 341     ostringstream err1;
 342     char reply;
 343     int nSeqs = alignPtr->getNumSeqs();
 344
 345     pathToRoot = new TreeNode*[nSeqs];
 346     distToNode = new float[nSeqs];
 347     dmat = new double*[nSeqs];
 348
 349     for (i = 0; i < nSeqs; i++)
 350     {
 351         dmat[i] = new double[nSeqs];
 352     }
 353
 354     if (nSeqs >= 2)
 355     {
 356         /*
 357          * for each leaf, determine all nodes between the leaf and the root;
 358          */
 359         for (i = 0; i < nSeqs; i++)
 360         {
 361             depth = 0;dist = 0.0;
 362             p = olptr[i];
 363             while (p != NULL)
 364             {
 365                 pathToRoot[depth] = p;
 366                 dist += p->dist;
 367                 distToNode[depth] = dist;
 368                 p = p->parent;
 369                 depth++;
 370             }
 371
 372             /*
 373              * for each pair....
 374              */
 375             for (j = 0; j < i; j++)
 376             {
 377                 p = olptr[j];
 378                 dist = 0.0;
 379                 /*
 380                  * find the common ancestor.
 381                  */
 382                 found = false;
 383                 n = 0;
 384                 while ((found == false) && (p->parent != NULL))
 385                 {
 386                     for (k = 0; k < depth; k++)
 387                     if (p->parent == pathToRoot[k])
 388                     {
 389                         found = true;
 390                         n = k;
 391                     }
 392
 393                     dist += p->dist;
 394                     p = p->parent;
 395                 }
 396
 397                 dmat[i][j] = dist + distToNode[n - 1];
 398             }
 399         }
 400
 401         nerrs = 0;
 402         for (i = 0; i < nSeqs; i++)
 403         {
 404             dmat[i][i] = 0.0;
 405             for (j = 0; j < i; j++)
 406             {
 407                 if (dmat[i][j] < 0.01)
 408                 {
 409                     dmat[i][j] = 0.01;
 410                 }
 411                 if (dmat[i][j] > 1.0)
 412                 {
 413                     if (dmat[i][j] > 1.1 && nerrs < MAXERRS)
 414                     {
 415                         seq1[nerrs] = i;
 416                         seq2[nerrs] = j;
 417                         badDist[nerrs] = dmat[i][j];
 418                         nerrs++;
 419                     }
 420                     dmat[i][j] = 1.0;
 421                 }
 422             }
 423         }
 424         if (nerrs > 0)
 425         {
 426             string errMess = "The following sequences are too divergent to be aligned:\n";
 427
 428             for (i = 0; i < nerrs && i < 5; i++)
 429             {
 430                 err1 << "           " << alignPtr->getName(seq1[i] + 1) << " and "
 431                      << alignPtr->getName(seq2[i] + 1) << " (distance "
 432                      << setprecision(3) << badDist[i] << ")\n";
 433             }
 434             errMess += err1.str();
 435             errMess += "(All distances should be between 0.0 and 1.0)\n";
 436             errMess += "This may not be fatal but you have been warned!\n";
 437             errMess += "SUGGESTION: Remove one or more problem sequences and try again";
 438
 439             if (clustalw::userParameters->getInteractive())
 440             {
 441                 reply = clustalw::utilityObject->promptForYesNo(errMess.c_str(), "Continue ");
 442             }
 443             else
 444             {
 445                 reply = 'y';
 446             }
 447             if ((reply != 'y') && (reply != 'Y'))
 448             {
 449                 return 0;
 450             }
 451         }
 452     }
 453     else
 454     {
 455         for (i = 0; i < nSeqs; i++)
 456         {
 457             for (j = 0; j < i; j++)
 458             {
 459                 dmat[i][j] = (*distMat)(i + 1, j + 1);
 460             }
 461         }
 462     }
 463     delete []pathToRoot;
 464     delete []distToNode;
 465
 466     double value;
 467     for (i = 0; i < nSeqs; i++)
 468     {
 469         distMat->SetAt(i + 1, i + 1, 0.0);
 470         for (j = 0; j < i; j++)
 471         {
 472             value = 100.0 - (dmat[i][j]) * 100.0;
 473             distMat->SetAt(i + 1, j + 1, value);
 474             distMat->SetAt(j + 1, i + 1, value);
 475         }
 476     }
 477
 478     for (i = 0; i < nSeqs; i++)
 479     {
 480         delete [] dmat[i];
 481     }
 482
 483     delete []dmat;
 484
 485     return 1;
 486
 487 }
 488
 489 /** *************************************************************************
 490  * Private functions!!!!!!!!!!!!!!!                                         *
 491  ****************************************************************************/
 492
 493 /**
 494  *
 495  * @param ptree
 496  * @param parent
 497  * @param file
 498  */
 499 void Tree::createTree(clustalw::TreeNode* ptree, clustalw::TreeNode* parent, ifstream* file)
 500 {
 501     TreeNode* p;
 502
 503     int i, type;
 504     float dist;
 505     string name;
 506
 507
 508     // is this a node or a leaf ?
 509     skipSpace(file);
 510     charFromFile = file->get();
 511     if (charFromFile == '(')
 512     {
 513         // this must be a node....
 514         type = NODE;
 515         name = "";
 516         ptrs[ntotal] = nptr[nnodes] = ptree;
 517         nnodes++;
 518         ntotal++;
 519
 520         createNode(ptree, parent);
 521
 522         p = ptree->left;
 523         createTree(p, ptree, file);
 524
 525         if (charFromFile == ',')
 526         {
 527             p = ptree->right;
 528             createTree(p, ptree, file);
 529             if (charFromFile == ',')
 530             {
 531                 ptree = insertNode(ptree);
 532                 ptrs[ntotal] = nptr[nnodes] = ptree;
 533                 nnodes++;
 534                 ntotal++;
 535                 p = ptree->right;
 536                 createTree(p, ptree, file);
 537                 rootedTree = false;
 538             }
 539         }
 540
 541         skipSpace(file);
 542         charFromFile = file->get();
 543     }
 544
 545     // ...otherwise, this is a leaf
 546
 547     else
 548     {
 549         type = LEAF;
 550         ptrs[ntotal++] = lptr[numSeq++] = ptree;
 551         // get the sequence name
 552         name = "";
 553         name += charFromFile;
 554         charFromFile = file->get();
 555
 556         i = 1;
 557         while ((charFromFile != ':') && (charFromFile != ',') && (charFromFile != ')'))
 558         {
 559             if (i < MAXNAMES)
 560             {
 561                 name += charFromFile;
 562                 i++;
 563             }
 564             charFromFile = file->get();
 565         }
 566
 567         if (charFromFile != ':')
 568         {
 569             clustalw::userParameters->setDistanceTree(false);
 570             dist = 0.0;
 571         }
 572     }
 573
 574     // get the distance information
 575
 576     dist = 0.0;
 577     if (charFromFile == ':')
 578     {
 579         skipSpace(file);
 580         (*file) >> dist;
 581         skipSpace(file);
 582         charFromFile = file->get();
 583     }
 584     setInfo(ptree, parent, type, name, dist);
 585 }
 586
 587 /**
 588  *
 589  * @param pptr
 590  * @param parent
 591  */
 592 void Tree::createNode(TreeNode* pptr, TreeNode* parent)
 593 {
 594     TreeNode* t;
 595     pptr->parent = parent;
 596     t = avail();
 597     pptr->left = t;
 598     t = avail();
 599     pptr->right = t;
 600 }
 601
 602 /**
 603  *
 604  * @param pptr
 605  * @return
 606  */
 607 TreeNode* Tree::insertNode(TreeNode* pptr)
 608 {
 609
 610     TreeNode* newnode;
 611
 612     newnode = avail();
 613     createNode(newnode, pptr->parent);
 614
 615     newnode->left = pptr;
 616     pptr->parent = newnode;
 617
 618     setInfo(newnode, pptr->parent, NODE, "", 0.0);
 619
 620     return newnode;
 621 }
 622
 623 /**
 624  *
 625  * @param p
 626  */
 627 void Tree::clearTree(clustalw::TreeNode* p)
 628 {
 629     clearTreeNodes(p);
 630     delete [] nptr;
 631     nptr  = NULL;
 632     delete [] ptrs;
 633     ptrs  = NULL;
 634     delete [] lptr;
 635     lptr  = NULL;
 636     delete [] olptr;
 637     olptr  = NULL;
 638 }
 639
 640 /**
 641  *
 642  * @param p
 643  */
 644 void Tree::clearTreeNodes(clustalw::TreeNode* p)
 645 {
 646     if (p == NULL)
 647     {
 648         p = root;
 649     }
 650     if (p->left != NULL)
 651     {
 652         clearTreeNodes(p->left);
 653     }
 654     if (p->right != NULL)
 655     {
 656         clearTreeNodes(p->right);
 657     }
 658     p->left = NULL;
 659     p->right = NULL;
 660
 661     delete p;
 662     p  = NULL;
 663 }
 664
 665
 666
 667
 668 /**
 669  *
 670  * @param
 671  * @param
 672  * @return
 673  */
 674 void Tree::debugPrintAllNodes(int nseqs)
 675 {
 676   clustalw::TreeNode *p;
 677   int i;
 678    float diff, maxDist;
 679
 680     cerr << "\nDEBUG: reportAllNodes\n";
 681     for (i = 0; i < ntotal; i++) {
 682         p = ptrs[i];
 683         //        ios::sync_with_stdio();
 684
 685         // same design as TreeNode
 686         if (p->parent == NULL)
 687             diff = calcRootMean(p, &maxDist);
 688         else
 689             diff = calcMean(p, &maxDist, nseqs);
 690         fprintf(stdout,
 691                 "i=%d p=%p: parent=%p left=%p right=%p dist=%f diff=%f\n",
 692                 i, (void*)p, (void*)p->parent, (void*)p->left, (void*)p->right,
 693                 p->dist, diff);
 694     }
 695 }
 696
 697
 698
 699
 700
 701 /**
 702  *
 703  * @param ptree
 704  * @param nseqs
 705  * @return
 706  */
 707 clustalw::TreeNode* Tree::reRoot(clustalw::TreeNode* ptree, int nseqs)
 708 {
 709     clustalw::TreeNode *p, *rootNode, *rootPtr;
 710     float diff, minDiff = 0.0, minDepth = 1.0, maxDist;
 711     int i;
 712     bool first = true;
 713
 714     // find the difference between the means of leaf->node
 715     // distances on the left and on the right of each node
 716     rootPtr = ptree;
 717     for (i = 0; i < ntotal; i++)
 718     {
 719         p = ptrs[i];
 720         if (p->parent == NULL)
 721         {
 722             /* AW Bug 94: p->parent must be chosen as rootNode
 723                (non-optimized executables (-O0) never do), otherwise
 724                insertRoot fails.
 725                Is p->parent == NULL valid at all?
 726                Simplest thing for now is to continue here. Tree code
 727                needs serious dismantling anyway. See debugPrintAllNodes
 728             */
 729
 730             continue;
 731             //diff = calcRootMean(p, &maxDist);
 732         }
 733         else
 734         {
 735             diff = calcMean(p, &maxDist, nseqs);
 736         }
 737
 738         if ((diff == 0) || ((diff > 0) && (diff < 2 *p->dist)))
 739         {
 740             if ((maxDist < minDepth) || (first == true))
 741             {
 742                 first = false;
 743                 rootPtr = p;
 744                 minDepth = maxDist;
 745                 minDiff = diff;
 746             }
 747         }
 748
 749     }
 750
 751
 752     // insert a new node as the ancestor of the node which produces the shallowest
 753     // tree.
 754     /* AW Bug 94: could also be prevented here */
 755     if (rootPtr == ptree)
 756     {
 757         minDiff = rootPtr->left->dist + rootPtr->right->dist;
 758         rootPtr = rootPtr->right;
 759     }
 760     rootNode = insertRoot(rootPtr, minDiff);
 761
 762     diff = calcRootMean(rootNode, &maxDist);
 763
 764     return rootNode;
 765 }
 766
 767 /**
 768  *
 769  * @param p
 770  * @param diff
 771  * @return
 772  */
 773 clustalw::TreeNode* Tree::insertRoot(clustalw::TreeNode* p, float diff)
 774 {
 775     clustalw::TreeNode *newp, *prev, *q, *t;
 776     float dist, prevDist, td;
 777
 778     newp = avail();
 779
 780     if (p->parent==NULL) {
 781         // AW bug 94: question remains if access here should be handled differently
 782         cerr << "\n\n*** INTERNAL ERROR: Tree::insertRoot: TreeNode p->parent is NULL\n";
 783         cerr << "To help us fix this bug, please send sequence file and used options to clustalw@ucd.ie\n";
 784         exit(1);
 785     }
 786
 787     t = p->parent;
 788     prevDist = t->dist;
 789
 790     p->parent = newp;
 791
 792     dist = p->dist;
 793
 794     p->dist = diff / 2;
 795
 796     if (p->dist < 0.0)
 797     {
 798         p->dist = 0.0;
 799     }
 800
 801     if (p->dist > dist)
 802     {
 803         p->dist = dist;
 804     }
 805
 806     t->dist = dist - p->dist;
 807
 808     newp->left = t;
 809     newp->right = p;
 810     newp->parent = NULL;
 811     newp->dist = 0.0;
 812     newp->leaf = NODE;
 813
 814     if (t->left == p)
 815     {
 816         t->left = t->parent;
 817     }
 818     else
 819     {
 820         t->right = t->parent;
 821     }
 822
 823     prev = t;
 824     q = t->parent;
 825
 826     t->parent = newp;
 827
 828     while (q != NULL)
 829     {
 830         if (q->left == prev)
 831         {
 832             q->left = q->parent;
 833             q->parent = prev;
 834             td = q->dist;
 835             q->dist = prevDist;
 836             prevDist = td;
 837             prev = q;
 838             q = q->left;
 839         }
 840         else
 841         {
 842             q->right = q->parent;
 843             q->parent = prev;
 844             td = q->dist;
 845             q->dist = prevDist;
 846             prevDist = td;
 847             prev = q;
 848             q = q->right;
 849         }
 850     }
 851
 852     /*
 853      * remove the old root node
 854      */
 855     q = prev;
 856     if (q->left == NULL)
 857     {
 858         dist = q->dist;
 859         q = q->right;
 860         q->dist += dist;
 861         q->parent = prev->parent;
 862
 863         if (prev->parent->left == prev)
 864         {
 865             prev->parent->left = q;
 866         }
 867         else
 868         {
 869             prev->parent->right = q;
 870         }
 871
 872         prev->right = NULL;
 873     }
 874     else
 875     {
 876         dist = q->dist;
 877         q = q->left;
 878         q->dist += dist;
 879         q->parent = prev->parent;
 880
 881         if (prev->parent->left == prev)
 882         {
 883             prev->parent->left = q;
 884         }
 885         else
 886         {
 887             prev->parent->right = q;
 888         }
 889
 890         prev->left = NULL;
 891     }
 892
 893     return (newp);
 894 }
 895
 896 /**
 897  *
 898  * @param root
 899  * @param maxDist
 900  * @return
 901  */
 902 float Tree::calcRootMean(clustalw::TreeNode* root, float *maxDist)
 903 {
 904     float dist, leftSum = 0.0, rightSum = 0.0, leftMean, rightMean, diff;
 905     clustalw::TreeNode* p;
 906     int i;
 907     int numLeft, numRight;
 908     int direction;
 909
 910     // for each leaf, determine whether the leaf is left or right of the root.
 911
 912     dist = (*maxDist) = 0;
 913     numLeft = numRight = 0;
 914     for (i = 0; i < numSeq; i++)
 915     {
 916         p = lptr[i];
 917         dist = 0.0;
 918         while (p->parent != root)
 919         {
 920             dist += p->dist;
 921             p = p->parent;
 922         }
 923
 924         if (p == root->left)
 925         {
 926             direction = LEFT;
 927         }
 928         else
 929         {
 930             direction = RIGHT;
 931         }
 932
 933         dist += p->dist;
 934
 935         if (direction == LEFT)
 936         {
 937             leftSum += dist;
 938             numLeft++;
 939         }
 940         else
 941         {
 942             rightSum += dist;
 943             numRight++;
 944         }
 945         if (dist > (*maxDist))
 946         {
 947             *maxDist = dist;
 948         }
 949     }
 950
 951     leftMean = leftSum / numLeft;
 952     rightMean = rightSum / numRight;
 953
 954     diff = leftMean - rightMean;
 955     return diff;
 956 }
 957
 958 /**
 959  *
 960  * @param nptr
 961  * @param maxDist
 962  * @param nSeqs
 963  * @return
 964  */
 965 float Tree::calcMean(clustalw::TreeNode* nptr, float *maxDist, int nSeqs)
 966 {
 967     float dist, leftSum = 0.0, rightSum = 0.0, leftMean, rightMean, diff;
 968     clustalw::TreeNode* p;
 969     clustalw::TreeNode** pathToRoot;
 970     float *distToNode;
 971     int depth = 0, i, j, n = 0;
 972     int numLeft, numRight;
 973     int direction;
 974     bool found;
 975
 976     pathToRoot = new clustalw::TreeNode*[nSeqs];
 977     distToNode = new float[nSeqs];
 978     // determine all nodes between the selected node and the root;
 979
 980     depth = 0;
 981     (*maxDist) = dist = 0.0;
 982     numLeft = numRight = 0;
 983     p = nptr;
 984     while (p != NULL)
 985     {
 986         pathToRoot[depth] = p;
 987         dist += p->dist;
 988         distToNode[depth] = dist;
 989         p = p->parent;
 990         depth++;
 991     }
 992
 993     // for each leaf, determine whether the leaf is left or right of the node.
 994     // (RIGHT = descendant, LEFT = not descendant)
 995
 996     for (i = 0; i < numSeq; i++)
 997     {
 998         p = lptr[i];
 999         if (p == nptr)
1000         {
1001             direction = RIGHT;
1002             dist = 0.0;
1003         }
1004         else
1005         {
1006             direction = LEFT;
1007             dist = 0.0;
1008
1009             // find the common ancestor.
1010
1011             found = false;
1012             n = 0;
1013             while ((found == false) && (p->parent != NULL))
1014             {
1015                 for (j = 0; j < depth; j++)
1016                 if (p->parent == pathToRoot[j])
1017                 {
1018                     found = true;
1019                     n = j;
1020                 }
1021
1022                 dist += p->dist;
1023                 p = p->parent;
1024             }
1025             if (p == nptr)
1026             {
1027                 direction = RIGHT;
1028             }
1029         }
1030
1031         if (direction == LEFT)
1032         {
1033             leftSum += dist;
1034             leftSum += distToNode[n - 1];
1035             numLeft++;
1036         }
1037         else
1038         {
1039             rightSum += dist;
1040             numRight++;
1041         }
1042
1043         if (dist > (*maxDist))
1044         {
1045             *maxDist = dist;
1046         }
1047     }
1048
1049     delete [] distToNode;
1050     distToNode = NULL;
1051     delete [] pathToRoot;
1052     pathToRoot = NULL;
1053
1054     leftMean = leftSum / numLeft;
1055     rightMean = rightSum / numRight;
1056
1057     diff = leftMean - rightMean;
1058     return (diff);
1059 }
1060
1061 /**
1062  *
1063  */
1064 void Tree::orderNodes()
1065 {
1066     int i;
1067     clustalw::TreeNode* p;
1068
1069     for (i = 0; i < numSeq; i++)
1070     {
1071         p = lptr[i];
1072         while (p != NULL)
1073         {
1074             p->order++;
1075             p = p->parent;
1076         }
1077     }
1078 }
1079
1080 /**
1081  *
1082  * @param leaf
1083  * @return
1084  */
1085 int Tree::calcWeight(int leaf)
1086 {
1087     clustalw::TreeNode* p;
1088     float weight = 0.0;
1089
1090     p = olptr[leaf];
1091     while (p->parent != NULL)
1092     {
1093         weight += p->dist / p->order;
1094         p = p->parent;
1095     }
1096
1097     weight *= 100.0;
1098
1099     return ((int)weight);
1100 }
1101
1102 /**
1103  * skipSpace is used to skip all the spaces at the begining of a file. The next read
1104  * will give a character other than a space.
1105  * @param file
1106  */
1107 void Tree::skipSpace(ifstream* file)
1108 {
1109     char c;
1110
1111     do
1112     {
1113         c = file->get();
1114     }
1115     while (isspace(c));
1116
1117     file->putback(c);
1118 }
1119
1120 /**
1121  *
1122  * @param p
1123  * @param nextGroups
1124  * @param nSeqs
1125  * @param stepsPtr
1126  */
1127 void Tree::groupSeqs(clustalw::TreeNode* p, int *nextGroups, int nSeqs, AlignmentSteps* stepsPtr)
1128 {
1129     int i;
1130     int *tmpGroups;
1131
1132     tmpGroups = new int[nSeqs + 1];
1133     for (i = 0; i < nSeqs; i++)
1134     {
1135         tmpGroups[i] = 0;
1136     }
1137
1138     if (p->left != NULL)
1139     {
1140         if (p->left->leaf == NODE)
1141         {
1142             groupSeqs(p->left, nextGroups, nSeqs, stepsPtr);
1143
1144             for (i = 0; i < nSeqs; i++)
1145                 if (nextGroups[i] != 0)
1146                 {
1147                     tmpGroups[i] = 1;
1148                 }
1149         }
1150         else
1151         {
1152             markGroup1(p->left, tmpGroups, nSeqs);
1153         }
1154
1155     }
1156
1157     if (p->right != NULL)
1158     {
1159         if (p->right->leaf == NODE)
1160         {
1161             groupSeqs(p->right, nextGroups, nSeqs, stepsPtr);
1162             for (i = 0; i < nSeqs; i++)
1163                 if (nextGroups[i] != 0)
1164                 {
1165                     tmpGroups[i] = 2;
1166                 }
1167         }
1168         else
1169         {
1170             markGroup2(p->right, tmpGroups, nSeqs);
1171         }
1172         stepsPtr->saveSet(nSeqs, tmpGroups);
1173     }
1174
1175     for (i = 0; i < nSeqs; i++)
1176     {
1177         nextGroups[i] = tmpGroups[i];
1178     }
1179
1180     delete [] tmpGroups;
1181     tmpGroups  = NULL;
1182 }
1183
1184 /**
1185  *
1186  * @param p
1187  * @param groups
1188  * @param n
1189  */
1190 void Tree::markGroup1(clustalw::TreeNode* p, int *groups, int n)
1191 {
1192     int i;
1193
1194     for (i = 0; i < n; i++)
1195     {
1196         if (olptr[i] == p)
1197         {
1198             groups[i] = 1;
1199         }
1200         else
1201         {
1202             groups[i] = 0;
1203         }
1204     }
1205 }
1206
1207 /**
1208  *
1209  * @param p
1210  * @param groups
1211  * @param n
1212  */
1213 void Tree::markGroup2(clustalw::TreeNode* p, int *groups, int n)
1214 {
1215     int i;
1216
1217     for (i = 0; i < n; i++)
1218     {
1219         if (olptr[i] == p)
1220         {
1221             groups[i] = 2;
1222         }
1223         else if (groups[i] != 0)
1224         {
1225             groups[i] = 1;
1226         }
1227     }
1228 }
1229
1230 /**
1231  *
1232  * @return
1233  */
1234 clustalw::TreeNode* Tree::avail()
1235 {
1236     clustalw::TreeNode* p;
1237     p = new TreeNode;
1238     p->left = NULL;
1239     p->right = NULL;
1240     p->parent = NULL;
1241     p->dist = 0.0;
1242     p->leaf = 0;
1243     p->order = 0;
1244     return (p);
1245 }
1246
1247 /**
1248  *
1249  * @param p
1250  * @param parent
1251  * @param pleaf
1252  * @param pname
1253  * @param pdist
1254  */
1255 void Tree::setInfo(TreeNode* p, TreeNode* parent, int pleaf, string pname, float
1256                      pdist)
1257 {
1258     p->parent = parent;
1259     p->leaf = pleaf;
1260     p->dist = pdist;
1261     p->order = 0;
1262     p->name = pname;
1263     if (p->leaf == true)
1264     {
1265         p->left = NULL;
1266         p->right = NULL;
1267     }
1268 }
1269
1270 }