binaries/src/clustalw/src/tree/UPGMA/UPGMAAlgorithm.cpp

   1 #ifdef HAVE_CONFIG_H
   2     #include "config.h"
   3 #endif
   4 #include "UPGMAAlgorithm.h"
   5 #include <cmath>
   6 #include <ctime>
   7 #include <iomanip>
   8 #include <sstream>
   9
  10 #include "../../general/SymMatrix.h"
  11 #include "../../general/debuglogObject.h"
  12 #include "../../general/clustalw.h"
  13 #include "upgmadata.h"
  14
  15 namespace clustalw
  16 {
  17
  18 UPGMAAlgorithm::UPGMAAlgorithm()
  19 : overwriteMatrix(false),
  20   numSeqs(0),
  21   verbose(false),
  22   orderNode1(0),
  23   orderNode2(0),
  24   orderNewNode(0)
  25 {}
  26
  27 auto_ptr<AlignmentSteps> UPGMAAlgorithm::generateTree(RootedGuideTree* phyTree,
  28                                           DistMatrix* distMat,
  29                                           SeqInfo* seqInfo, bool overwrite,
  30                                           ofstream* tree)
  31 {
  32     if (tree == 0 || !tree->is_open())
  33     {
  34         verbose = false;
  35     }
  36
  37     if (verbose)
  38     {
  39         (*tree) << "\n\n\t\t\tUPGMA Method\n"
  40                 << "\n\n This is a ROOTED tree\n"
  41                 << "\n Numbers in parentheses are branch lengths\n\n";
  42     }
  43
  44     progSteps.reset(new AlignmentSteps);
  45
  46     Node** clusters;
  47     Node* root;
  48     numSeqs = seqInfo->numSeqs;
  49     const int sizeDistMat = ((numSeqs + 1) * (numSeqs + 2)) / 2;
  50
  51     double* elements = overwrite ?
  52                       distMat->getDistMatrix(seqInfo->firstSeq, seqInfo->numSeqs) :
  53                       (double *)memcpy(new double[sizeDistMat],
  54                       distMat->getDistMatrix(seqInfo->firstSeq, seqInfo->numSeqs),
  55                       sizeDistMat * sizeof(double));
  56
  57     clusters = initialiseNodes(elements, seqInfo->firstSeq);
  58     root = doUPGMA(clusters, tree);
  59
  60     phyTree->setRoot(root);
  61     delete [] clusters;
  62
  63     if(!overwrite)
  64     {
  65         delete [] elements;
  66     }
  67     distMat->clearSubArray();
  68
  69     return progSteps;
  70 }
  71
  72
  73 Node** UPGMAAlgorithm::initialiseNodes(double* distanceMatrix, int fSeq)
  74 {
  75     int firstSeq = fSeq;
  76
  77
  78     Node** nodes = new Node*[numSeqs];
  79     Node** nodeIter = nodes;
  80
  81     *nodes = new Node(firstSeq, 0, 0);
  82
  83     distanceMatrix++;
  84
  85     // Move to first position in distanceMatrix.
  86     for(int elementIndex = 1, e = numSeqs; elementIndex < e;
  87         distanceMatrix += ++elementIndex)
  88     {
  89         Node* newcluster = new Node(elementIndex + firstSeq,
  90                                     distanceMatrix, elementIndex);
  91         (*nodeIter++)->next = newcluster;
  92         *nodeIter = newcluster;
  93     }
  94
  95   return nodes;
  96 }
  97
  98 void UPGMAAlgorithm::printAllNodes(Node** nodes)
  99 {
 100     int numNodes = 0;
 101     for(Node* nodeIter = *nodes; nodeIter; nodeIter = nodeIter->next)
 102     {
 103         numNodes++;
 104         cout << "Node " << numNodes << "\n";
 105         nodeIter->printElements();
 106         cout << "\n\n";
 107     }
 108     cout << "There are " << numNodes << " nodes\n";
 109 }
 110
 111 Node* UPGMAAlgorithm::doUPGMA(Node** nodes, ofstream* tree)
 112 {
 113     if (tree == 0 || !tree->is_open())
 114     {
 115         verbose = false;
 116     }
 117
 118     string type1, type2;
 119     int step = 0;
 120
 121     while((*nodes)->next) // While there is more than 1 node.
 122     {
 123         step++;
 124         if (verbose)
 125         {
 126             (*tree) <<  "\n Cycle" << setw(4) << step << "     = ";
 127         }
 128         Node** ptrNodeWithMin;
 129
 130         /**
 131          * STEP 1 find node with min distance
 132          */
 133         ptrNodeWithMin = getNodeWithMinDist(nodes);
 134         Node* nodeToJoin2 = *ptrNodeWithMin;
 135         const int indexToMinDist = nodeToJoin2->indexToMinDist;
 136         Node* const nodeToJoin1 = nodes[indexToMinDist];
 137
 138         orderNode1 = nodeToJoin1->size;
 139         orderNode2 = nodeToJoin2->size;
 140         orderNewNode = orderNode1 + orderNode2;
 141
 142         /**
 143          * STEP 2 Recompute the dist Matrix row for nodeToJoin1
 144          * example: Join nodes 2 and 6
 145          *
 146          * 1  0
 147          * 2  X 0
 148          * 3  0 0 0
 149          * 4  0 0 0 0
 150          * 5  0 0 0 0 0
 151          * 6  0 0 0 0 0 0
 152          * 7  0 0 0 0 0 0 0
 153          */
 154         double* nodeToJoin2DistIter = nodeToJoin2->ptrToDistMatRow;
 155
 156         if (indexToMinDist != 0)
 157         {
 158             recomputeNodeToJoin1DistMatRow(nodeToJoin1, &nodeToJoin2DistIter);
 159         }
 160
 161         /**
 162          * STEP 3 Recompute all distances in column
 163          * example: Join nodes 2 and 6
 164          *
 165          * 1  0
 166          * 2  C 0
 167          * 3  0 X 0
 168          * 4  0 X 0 0
 169          * 5  0 X 0 0 0
 170          * 6  0 0 0 0 0 0
 171          * 7  0 X 0 0 0 0 0
 172          */
 173
 174         computeAllOtherDistsToNewNode(nodeToJoin1, nodeToJoin2, &nodeToJoin2DistIter);
 175
 176         /**
 177          * STEP 4 Add the step to the progSteps.
 178          * This creates the multiple alignment steps.
 179          */
 180         addAlignmentStep(&nodeToJoin1->allElements, &nodeToJoin2->allElements);
 181
 182
 183         double minDistance = (*ptrNodeWithMin)->minDist;
 184
 185         double height = 0.0;
 186         height = minDistance / 2.0;
 187
 188         if(verbose)
 189         {
 190             if(nodeToJoin1->allElements.size() > 1)
 191             {
 192                 type1 = "NODE: ";
 193             }
 194             else
 195             {
 196                 type1 = "SEQ: ";
 197             }
 198             if(nodeToJoin2->allElements.size() > 1)
 199             {
 200                 type2 = "NODE: ";
 201             }
 202             else
 203             {
 204                 type2 = "SEQ: ";
 205             }
 206             (*tree) << type1 << nodeToJoin1->allElements[0] << " (" << setw(9)
 207                     << setprecision(5) << height << ") joins " << type2
 208                     << setw(4) << nodeToJoin2->allElements[0] << " ("
 209                     << setw(9) << setprecision(5) << height << ")";
 210         }
 211         /**
 212          * STEP 5 merge 2 nodes
 213          */
 214         nodeToJoin1->merge(ptrNodeWithMin, height);
 215
 216     }
 217
 218     return *nodes;
 219 }
 220
 221 /**
 222  * This function returns a Node object that has the minimum distance of
 223  * all the nodes left.
 224  */
 225 Node** UPGMAAlgorithm::getNodeWithMinDist(Node** nodes)
 226 {
 227     Node** ptrMinNode = NULL;
 228     double minDistance = numeric_limits<double>::max();
 229     //Start at node 1, check see if it points to something, then move on the next one.
 230     Node** nodeIter;
 231     for(nodeIter = &((*nodes)->next); *nodeIter;
 232         nodeIter = &(*nodeIter)->next)
 233     {
 234         if ((*nodeIter)->getMinDist() < minDistance)
 235         {
 236             minDistance = (*nodeIter)->getMinDist();
 237             ptrMinNode = nodeIter; // ptrMinNode will hold a ptr to node with sm'st dist
 238         }
 239     }
 240     return ptrMinNode;
 241 }
 242
 243 /**
 244  * This function is used to recompute the nodeToJoin1 dist mat row. Each row in the distance
 245  * matrix has the distances to all nodes before it, not after. For example the dist mat row
 246  * for Node 3 would have the distances to nodes 1 and 2.
 247          * 1  0
 248          * 2  0 0
 249          * 3  X X 0      Dists to node 1 and 2.
 250          * 4  0 0 0 0
 251          * 5  0 0 0 0 0
 252          * 6  0 0 0 0 0 0
 253          * 7  0 0 0 0 0 0 0
 254  */
 255 void UPGMAAlgorithm::recomputeNodeToJoin1DistMatRow(Node* nodeToJoin1,
 256                                                     double** nodeToJoin2DistIter)
 257 {
 258     double* nodeToJoin1DistIter = nodeToJoin1->ptrToDistMatRow;
 259     // calculate new distance
 260     *nodeToJoin1DistIter = calcNewDist(*nodeToJoin1DistIter, **nodeToJoin2DistIter);
 261
 262     const double* minIndex1 = nodeToJoin1DistIter;
 263     nodeToJoin1DistIter++;
 264     (*nodeToJoin2DistIter)++;
 265     int numDistToUpdate = nodeToJoin1->numDists - 1;
 266
 267     // For each of the distances in nodeToJoin1
 268     while(numDistToUpdate > 0)
 269     {
 270         if (*nodeToJoin1DistIter >= 0)
 271         {
 272             // Calculate the average
 273             *nodeToJoin1DistIter = calcNewDist(*nodeToJoin1DistIter, **nodeToJoin2DistIter);
 274
 275             if (*nodeToJoin1DistIter < *minIndex1)
 276             {
 277                 minIndex1 = nodeToJoin1DistIter;
 278             }
 279         }
 280         nodeToJoin1DistIter++;
 281         (*nodeToJoin2DistIter)++;
 282         numDistToUpdate--;
 283     }
 284
 285     // We have found the minimum distance
 286     nodeToJoin1->minDist = *minIndex1;
 287     nodeToJoin1->indexToMinDist = minIndex1 - nodeToJoin1->ptrToDistMatRow;
 288 }
 289
 290 /**
 291  * This function is used to recompute all the other distances. It does this by calling
 292  * two other functions. We first compute all the distances until we get to node 2.
 293  * Then we call another function to do the rest. This is because nodes after node 2 may
 294  * have had EITHER node1 of node2 as their min distance.
 295  */
 296 void UPGMAAlgorithm::computeAllOtherDistsToNewNode(Node* nodeToJoin1, Node* nodeToJoin2,
 297                                                    double** nodeToJoin2DistIter)
 298 {
 299     computeDistsUpToNodeToJoin2(nodeToJoin1, nodeToJoin2, nodeToJoin2DistIter);
 300     computeDistsForNodesAfterNode2(nodeToJoin2);
 301 }
 302
 303 /**
 304  * STEP 3A:
 305  * This function recomputes the distances in column until we get to nodeToJoin2
 306  * example: Join nodes 2 and 6
 307  *
 308  * 1  0
 309  * 2  C 0
 310  * 3  0 X 0
 311  * 4  0 X 0 0
 312  * 5  0 X 0 0 0
 313  * 6  0 0 0 0 0 0
 314  * 7  0 0 0 0 0 0 0
 315  *
 316  * For each node until we get to nodeToJoin2
 317  *    If newdistance is less than the old min distance, set it to this.
 318  *    else if its greater than old minDist and indexTominDist is the same, recompute min
 319  *    else leave the min the same as it hasnt changed.
 320  */
 321 void UPGMAAlgorithm::computeDistsUpToNodeToJoin2(Node* nodeToJoin1, Node* nodeToJoin2, double** nodeToJoin2DistIter)
 322 {
 323     const int indexToMinDist = nodeToJoin2->indexToMinDist;
 324
 325     movePtrPastUnusedDistances(nodeToJoin2DistIter);
 326
 327     Node* nodeIter;
 328     // For each node until we get to the second node we are joining
 329     for(nodeIter = nodeToJoin1->next; nodeIter != nodeToJoin2; nodeIter = nodeIter->next)
 330     {
 331         (*nodeToJoin2DistIter)++; // Skip the distance to the node we are joining with
 332         movePtrPastUnusedDistances(nodeToJoin2DistIter);
 333
 334         double distToNode = nodeIter->ptrToDistMatRow[indexToMinDist];
 335
 336         double newDistToNode = calcNewDist(distToNode, **nodeToJoin2DistIter);
 337         nodeIter->ptrToDistMatRow[indexToMinDist] = newDistToNode;
 338
 339
 340         if (newDistToNode < nodeIter->minDist)
 341         { /** new value is smaller than current min. */
 342             nodeIter->minDist = newDistToNode;
 343             nodeIter->indexToMinDist = indexToMinDist;
 344         }
 345         else if ((newDistToNode > nodeIter->minDist) &&
 346                  (nodeIter->indexToMinDist == indexToMinDist))
 347         { /** indexToMinDist was the min dist, but it is now a larger num, recompute */
 348             nodeIter->findMinDist();
 349         }
 350
 351     }
 352
 353 }
 354
 355 /**
 356  * STEP 3B Recompute distance for nodes after nodeToJoin2
 357  * example: Join nodes 2 and 6
 358  *
 359  * 1  0
 360  * 2  C 0
 361  * 3  0 C 0
 362  * 4  0 C 0 0
 363  * 5  0 C 0 0 0
 364  * 6  0 0 0 0 0 0
 365  * 7  0 X 0 0 0 0 0
 366  *
 367  * For each node until we get to the end.
 368  *     if dist is less than minDist, set mindist to new distance, set index to index
 369  *     else if dist is greater than mindist and the index was either nodetojoin1 or
 370  *     nodetojoin2, recompute distance, set entry for nodetojoin2 to NOTUSED
 371  *     else set the distance to unused.
 372  */
 373
 374 void UPGMAAlgorithm::computeDistsForNodesAfterNode2(Node* nodeToJoin2)
 375 {
 376     int indexToNode2 = nodeToJoin2->numDists;
 377     const int indexToMinDist = nodeToJoin2->indexToMinDist;
 378
 379     Node* nodeIter;
 380
 381     for(nodeIter = nodeToJoin2->next; nodeIter; nodeIter = nodeIter->next)
 382     {
 383         double &distUpdate = nodeIter->ptrToDistMatRow[indexToMinDist];
 384
 385         distUpdate =
 386           ((distUpdate * orderNode1) +
 387            (nodeIter->ptrToDistMatRow[indexToNode2] * orderNode2))
 388           / orderNewNode;
 389
 390         /* The comparison (distUpdate > nodeIter->minDist) is unsafe.
 391          * Specifically, we get different results for 32/64bit machines,
 392          * which leads to different branching in the if/else statement
 393          * and nasty behaviour down the line.
 394          * Using all of Pfam as a benchmark distUpdate can 'wobble' by 40 ULPs
 395          * (Unit of Least Precision) which is difficult to translate into
 396          * a maximum relative error, so we pick COMPARE_REL_EPSILON
 397          * phenomenologically: approx 2E-15 was the biggest we saw in Pfam,
 398          * but suggest 1E-14 (for good measure).
 399          * Using ULPs, eg, *(long long int*)&(distUpdate),
 400          * would be better (more elegant?) but gives problems
 401          * with aliasing and/or performance reduction.
 402          * FS, 2009-04-30
 403          */
 404 #define COMPARE_REL_EPSILON 1E-14
 405         if ( (distUpdate < nodeIter->minDist) &&
 406              ((nodeIter->minDist-distUpdate) /
 407               nodeIter->minDist > COMPARE_REL_EPSILON) )
 408           { /** new distance is smaller */
 409             nodeIter->minDist = distUpdate;
 410             nodeIter->indexToMinDist = indexToMinDist;
 411           }
 412         else if (((distUpdate > nodeIter->minDist) &&
 413                   ((distUpdate-nodeIter->minDist) /
 414                    distUpdate > COMPARE_REL_EPSILON) &&
 415                   (nodeIter->indexToMinDist == indexToMinDist)) ||
 416                  (nodeIter->indexToMinDist == indexToNode2))
 417           { /** if old min dist was to either nodeToJoin1 or nodeToJoin2 */
 418             nodeIter->ptrToDistMatRow[indexToNode2] = BLANKDIST;
 419             nodeIter->findMinDist();
 420           }
 421         else
 422           {
 423             nodeIter->ptrToDistMatRow[indexToNode2] = BLANKDIST;
 424           }
 425     }
 426
 427 }
 428
 429 void UPGMAAlgorithm::addAlignmentStep(vector<int>* group1, vector<int>* group2)
 430 {
 431     int sizeGroup1 = group1->size();
 432     int sizeGroup2 = group2->size();
 433
 434     vector<int> groups;
 435     groups.resize(numSeqs + 1, 0);
 436     int sizeGroup = groups.size();
 437
 438     for(int i = 0; i < sizeGroup1 && (*group1)[i] < sizeGroup; i++)
 439     {
 440         groups[(*group1)[i]] = 1; // Set to be apart of group1
 441     }
 442
 443     for(int i = 0; i < sizeGroup2 && (*group2)[i] < sizeGroup; i++)
 444     {
 445         groups[(*group2)[i]] = 2; // Set to be apart of group1
 446     }
 447
 448     progSteps->saveSet(&groups);
 449 }
 450
 451 /**
 452  * At the moment we are only using average distance, UPGMA, but we can also use this
 453  * function to have different distance measures, min, max etc.
 454  */
 455 double UPGMAAlgorithm::calcNewDist(double dist1, double dist2)
 456 {
 457     return ((dist1 * orderNode1) + (dist2 * orderNode2)) / orderNewNode;
 458 }
 459
 460 }