binaries/src/clustalw/src/multipleAlign/MSA.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 /**
   7  * @author Mark Larkin, Conway Institute, UCD. mark.larkin@ucd.ie
   8  * Changes:
   9  *
  10  * Mark: 23-01-2007: There was a problem with running an alignment with only
  11  * one sequence in it. I needed to make a change to the multiSeqAlign function.
  12  ****************************************************************************/
  13 #ifdef HAVE_CONFIG_H
  14     #include "config.h"
  15 #endif
  16 #include "MSA.h"
  17 #include "MyersMillerProfileAlign.h"
  18 //#include "../phylogeneticTree/ClusterTree.h"
  19 #include "../general/debuglogObject.h"
  20
  21 namespace clustalw
  22 {
  23
  24 /**
  25  *
  26  * @param alnPtr
  27  * @param distMat
  28  * @param iStart
  29  * @param phylipName
  30  * @return
  31  */
  32 int MSA::multiSeqAlign(Alignment* alnPtr, DistMatrix* distMat, vector<int>* seqWeight, AlignmentSteps* progSteps, int iStart)
  33 {
  34
  35     if(!progSteps)
  36     {
  37         return 0;
  38     }
  39
  40     int* aligned;
  41     vector<int> group;
  42     int ix;
  43
  44     int* maxid;
  45     int max = 0, sum = 0;
  46     vector<int> treeWeight;
  47     int i = 0, j = 0, set = 0, iseq = 0;
  48     int entries = 0;
  49     int score = 0;
  50     int _numSteps = 0;
  51     utilityObject->info("Start of Multiple Alignment\n");
  52
  53     int _numSeqs = alnPtr->getNumSeqs();
  54
  55     vector<int> newOutputIndex(_numSeqs);
  56
  57     alnPtr->addSeqWeight(seqWeight);
  58
  59     ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
  60     _numSteps = progSteps->getNumSteps();
  61     // for each sequence, find the most closely related sequence
  62
  63     maxid = new int[_numSeqs + 1];
  64
  65     for (i = 1; i <= _numSeqs; i++)
  66     {
  67         maxid[i] =  -1;
  68         for (j = 1; j <= _numSeqs; j++)
  69             if (j != i && maxid[i] < (*distMat)(i, j))
  70             {
  71                 maxid[i] = static_cast<int>((*distMat)(i, j));
  72             }
  73     }
  74
  75     // group the sequences according to their relative divergence
  76
  77     if (iStart == 0)
  78     {
  79
  80         // start the multiple alignments.........
  81
  82         utilityObject->info("Aligning...");
  83         // first pass, align closely related sequences first....
  84
  85
  86         ix = 0;
  87         aligned = new int[_numSeqs + 1];
  88
  89         for (i = 0; i <= _numSeqs; i++)
  90         {
  91             aligned[i] = 0;
  92         }
  93
  94         const vector<vector<int> >* ptrToSets = progSteps->getSteps();
  95
  96
  97         for (set = 1; set <= _numSteps; ++set)
  98         {
  99             entries = 0;
 100             for (i = 1; i <= _numSeqs; i++)
 101             {
 102                 if (((*ptrToSets)[set][i] != 0) &&
 103                     (maxid[i] > userParameters->getDivergenceCutoff()))
 104                 {
 105                     entries++;
 106                     if (aligned[i] == 0)
 107                     {
 108                         if (userParameters->getOutputOrder() == INPUT)
 109                         {
 110                             ++ix;
 111                             newOutputIndex[i - 1] = i;
 112                         }
 113                         else
 114                         {
 115                             if(ix >= (int)newOutputIndex.size())
 116                             {
 117                                 cerr << "ERROR: size = " << newOutputIndex.size()
 118                                      << "ix = " << ix << "\n";
 119                                 exit(1);
 120                             }
 121                             else
 122                             {
 123                                 newOutputIndex[ix] = i;
 124                                 ++ix;
 125                             }
 126                         }
 127
 128                         aligned[i] = 1;
 129                     }
 130                 }
 131             }
 132
 133             if (entries > 0)
 134             {
 135                 #if DEBUGFULL
 136                     if(logObject && DEBUGLOG)
 137                     {
 138                         logObject->logMsg("Doing profile align");
 139                     }
 140                 #endif
 141                 score = alignAlgorithm->profileAlign(alnPtr, distMat, progSteps->getStep(set),
 142                                                      aligned);
 143             }
 144             else
 145             {
 146                 score = 0;
 147             }
 148
 149
 150             // negative score means fatal error... exit now!
 151
 152             if (score < 0)
 153             {
 154                 return (-1);
 155             }
 156             if(userParameters->getDisplayInfo())
 157             {
 158                 if ((entries > 0) && (score > 0))
 159                 {
 160                     utilityObject->info("Group %d: Sequences:%4d      Score:%d",
 161                                         set, entries, score);
 162                 }
 163                 else
 164                 {
 165                     utilityObject->info("Group %d:                     Delayed", set);
 166                 }
 167             }
 168         }
 169     }
 170
 171     else
 172     {
 173         aligned = new int[_numSeqs + 1];
 174         ix = 0;
 175         for (i = 1; i <= iStart + 1; i++)
 176         {
 177             aligned[i] = 1;
 178             ++ix;
 179             newOutputIndex[i - 1] = i;
 180         }
 181         for (i = iStart + 2; i <= _numSeqs; i++)
 182         {
 183             aligned[i] = 0;
 184         }
 185     }
 186
 187     // second pass - align remaining, more divergent sequences.....
 188
 189     // if not all sequences were aligned, for each unaligned sequence,
 190     // find it's closest pair amongst the aligned sequences.
 191
 192     group.resize(_numSeqs + 1);
 193     treeWeight.resize(_numSeqs);
 194
 195     for (i = 0; i < _numSeqs; i++)
 196     {
 197         treeWeight[i] = (*seqWeight)[i];
 198     }
 199
 200     // if we haven't aligned any sequences, in the first pass - align the
 201     // two most closely related sequences now
 202     if (ix == 0)
 203     {
 204         max =  -1;
 205         iseq = 0;
 206
 207         for (i = 1; i <= _numSeqs; i++)
 208         {
 209             for (j = i + 1; j <= _numSeqs; j++)
 210             {
 211                 if (max < (*distMat)(i, j))
 212                 {
 213                     max = static_cast<int>((*distMat)(i, j)); // Mark change 17-5-07
 214                     iseq = i;
 215                 }
 216             }
 217         }
 218         aligned[iseq] = 1;
 219         if (userParameters->getOutputOrder() == INPUT)
 220         {
 221             ++ix;
 222             newOutputIndex[iseq - 1] = iseq;
 223         }
 224         else
 225         {
 226             newOutputIndex[ix] = iseq;
 227             ++ix;
 228         }
 229     }
 230
 231     while (ix < _numSeqs)
 232     {
 233         for (i = 1; i <= _numSeqs; i++)
 234         {
 235             if (aligned[i] == 0)
 236             {
 237                 maxid[i] =  - 1;
 238                 for (j = 1; j <= _numSeqs; j++)
 239                     if ((maxid[i] < (*distMat)(i, j)) && (aligned[j] != 0))
 240                     {
 241                         maxid[i] = static_cast<int>((*distMat)(i, j));// Mark change 17-5-07
 242                     }
 243
 244             }
 245         }
 246         // find the most closely related sequence to those already aligned
 247
 248         max =  - 1;
 249         iseq = 0;
 250         for (i = 1; i <= _numSeqs; i++)
 251         {
 252             if ((aligned[i] == 0) && (maxid[i] > max))
 253             {
 254                 max = maxid[i];
 255                 iseq = i;
 256             }
 257         }
 258
 259
 260         // align this sequence to the existing alignment
 261         // weight sequences with percent identity with profile
 262         // OR...., multiply sequence weights from tree by percent identity with new sequence
 263         if (userParameters->getNoWeights() == false)
 264         {
 265             for (j = 0; j < _numSeqs; j++)
 266                 if (aligned[j + 1] != 0)
 267                 {
 268                     (*seqWeight)[j] = static_cast<int>(treeWeight[j] * (*distMat)(j + 1, iseq));
 269                 }
 270
 271             // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
 272
 273             sum = 0;
 274             for (j = 0; j < _numSeqs; j++)
 275             {
 276                 if (aligned[j + 1] != 0)
 277                 {
 278                     sum += (*seqWeight)[j];
 279                 }
 280             }
 281
 282             if (sum == 0)
 283             {
 284                 for (j = 0; j < _numSeqs; j++)
 285                 {
 286                     (*seqWeight)[j] = 1;
 287                 }
 288
 289                 sum = j;
 290             }
 291             for (j = 0; j < _numSeqs; j++)
 292             {
 293                 if (aligned[j + 1] != 0)
 294                 {
 295                     (*seqWeight)[j] = ((*seqWeight)[j] * INT_SCALE_FACTOR) / sum;
 296                     if ((*seqWeight)[j] < 1)
 297                     {
 298                         (*seqWeight)[j] = 1;
 299                     }
 300                 }
 301             }
 302         }
 303
 304         entries = 0;
 305         for (j = 1; j <= _numSeqs; j++)
 306         {
 307             if (aligned[j] != 0)
 308             {
 309                 group[j] = 1;
 310                 entries++;
 311             }
 312             else if (iseq == j)
 313             {
 314                 group[j] = 2;
 315                 entries++;
 316             }
 317         }
 318
 319         alnPtr->addSeqWeight(seqWeight);
 320         aligned[iseq] = 1;
 321
 322         score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
 323
 324         if (userParameters->getOutputOrder() == INPUT)
 325         {
 326             ++ix;
 327             newOutputIndex[iseq - 1] = iseq;
 328         }
 329         else
 330         {
 331             newOutputIndex[ix] = iseq;
 332             ++ix;
 333         }
 334     }
 335
 336     alnPtr->addOutputIndex(&newOutputIndex);
 337
 338     if(userParameters->getDisplayInfo())
 339     {
 340       int alignmentScore = alnPtr->alignScore(); // ?? check, FS, 2009-05-18
 341     }
 342
 343     delete alignAlgorithm;
 344     delete [] aligned;
 345     delete [] maxid;
 346     return (_numSeqs);
 347 }
 348
 349
 350 /**
 351  *
 352  * @param alnPtr
 353  * @param distMat
 354  * @param iStart
 355  * @param phylipName
 356  * @return
 357  */
 358 int MSA::seqsAlignToProfile(Alignment* alnPtr, DistMatrix* distMat, vector<int>* seqWeight, int iStart,
 359                            string phylipName)
 360 {
 361     int *aligned;
 362     vector<int> treeWeight;
 363     vector<int> group;
 364     int ix;
 365
 366     int *maxid;
 367     int max = 0;
 368     int i = 0, j = 0, iseq = 0;
 369     int sum = 0, entries = 0;
 370     int score = 0;
 371     int _numSeqs = alnPtr->getNumSeqs();
 372
 373     utilityObject->info("Start of Multiple Alignment\n");
 374
 375     ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
 376
 377     // calculate sequence weights according to branch lengths of the tree -
 378      // weights in global variable seq_weight normalised to sum to 100
 379     vector<int> newOutputIndex(_numSeqs);
 380     //groupTree.calcSeqWeights(0, _numSeqs, &seqWeight);
 381
 382     treeWeight.resize(_numSeqs);
 383     for (i = 0; i < _numSeqs; i++)
 384     {
 385         treeWeight[i] = (*seqWeight)[i];
 386     }
 387
 388     // for each sequence, find the most closely related sequence
 389
 390     maxid = new int[_numSeqs + 1];
 391
 392     for (i = 1; i <= _numSeqs; i++)
 393     {
 394         maxid[i] =  - 1;
 395         for (j = 1; j <= _numSeqs; j++)
 396         {
 397             if (maxid[i] < (*distMat)(i, j))
 398             {
 399                 maxid[i] = static_cast<int>((*distMat)(i, j)); // Mark change 17-5-07
 400             }
 401         }
 402     }
 403
 404     aligned = new int[_numSeqs + 1];
 405     ix = 0;
 406     for (i = 1; i <= iStart + 1; i++)
 407     {
 408         aligned[i] = 1;
 409         ++ix;
 410         newOutputIndex[i - 1] = i;
 411     }
 412     for (i = iStart + 2; i <= _numSeqs; i++)
 413     {
 414         aligned[i] = 0;
 415     }
 416
 417     // for each unaligned sequence, find it's closest pair amongst the
 418     // aligned sequences.
 419
 420     group.resize(_numSeqs + 1);
 421
 422     while (ix < _numSeqs)
 423     {
 424         if (ix > 0)
 425         {
 426             for (i = 1; i <= _numSeqs; i++)
 427             {
 428                 if (aligned[i] == 0)
 429                 {
 430                     maxid[i] =  - 1;
 431                     for (j = 1; j <= _numSeqs; j++)
 432                     {
 433                         if ((maxid[i] < (*distMat)(i, j)) && (aligned[j] != 0))
 434                         {
 435                             maxid[i] = static_cast<int>((*distMat)(i, j));
 436                         }
 437                     }
 438                 }
 439             }
 440         }
 441
 442         // find the most closely related sequence to those already aligned
 443
 444         max =  -1;
 445         for (i = 1; i <= _numSeqs; i++)
 446         {
 447             if ((aligned[i] == 0) && (maxid[i] > max))
 448             {
 449                 max = maxid[i];
 450                 iseq = i;
 451             }
 452         }
 453
 454         // align this sequence to the existing alignment
 455         entries = 0;
 456         for (j = 1; j <= _numSeqs; j++)
 457         {
 458             if (aligned[j] != 0)
 459             {
 460                 group[j] = 1;
 461                 entries++;
 462             }
 463             else if (iseq == j)
 464             {
 465                 group[j] = 2;
 466                 entries++;
 467             }
 468         }
 469
 470         aligned[iseq] = 1;
 471
 472         // multiply sequence weights from tree by percent
 473         // identity with new sequence
 474
 475         for (j = 0; j < _numSeqs; j++)
 476         {
 477             (*seqWeight)[j] = static_cast<int>(treeWeight[j] * (*distMat)(j + 1, iseq));
 478         }
 479
 480         //
 481         // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
 482         //
 483
 484         sum = 0;
 485         for (j = 0; j < _numSeqs; j++)
 486         {
 487             if (group[j + 1] == 1)
 488             {
 489                 sum += (*seqWeight)[j];
 490             }
 491         }
 492
 493         if (sum == 0)
 494         {
 495             for (j = 0; j < _numSeqs; j++)
 496             {
 497                 (*seqWeight)[j] = 1;
 498             }
 499
 500             sum = j;
 501         }
 502         for (j = 0; j < _numSeqs; j++)
 503         {
 504             (*seqWeight)[j] = ((*seqWeight)[j] * INT_SCALE_FACTOR) / sum;
 505             if ((*seqWeight)[j] < 1)
 506             {
 507                 (*seqWeight)[j] = 1;
 508             }
 509         }
 510         // Add the seqWeights to the Alignment object!!!!
 511         alnPtr->addSeqWeight(seqWeight);
 512
 513         score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
 514         utilityObject->info("Sequence:%d     Score:%d", iseq, score);
 515         if (userParameters->getOutputOrder() == INPUT)
 516         {
 517             ++ix;
 518             newOutputIndex[iseq - 1] = iseq;
 519         }
 520         else
 521         {
 522             newOutputIndex[ix] = iseq;
 523             ++ix;
 524         }
 525     }
 526
 527     delete [] aligned;
 528     delete [] maxid;
 529     delete alignAlgorithm;
 530     alnPtr->addOutputIndex(&newOutputIndex);
 531
 532     if(userParameters->getDisplayInfo())
 533     {
 534         alnPtr->alignScore();
 535     }
 536
 537     return (_numSeqs);
 538 }
 539
 540
 541 /**
 542  *
 543  * @param alnPtr
 544  * @param distMat
 545  * @return
 546  */
 547 int MSA::calcPairwiseForProfileAlign(Alignment* alnPtr, DistMatrix* distMat)
 548 {
 549     //Tree groupTree;
 550     int i, j, temp;
 551     int entries;
 552     int* aligned;
 553     vector<int> group;
 554     vector<int> seqWeight;
 555     float dscore;
 556     int score;
 557     int _numSeqs = alnPtr->getNumSeqs();
 558
 559     seqWeight.resize(_numSeqs);
 560     ProfileAlignAlgorithm* alignAlg = new MyersMillerProfileAlign;
 561
 562     utilityObject->info("Start of Initial Alignment");
 563     /* calculate sequence weights according to branch lengths of the tree -
 564      * weights in global variable seq_weight normalised to sum to INT_SCALE_FACTOR */
 565
 566     temp = INT_SCALE_FACTOR / _numSeqs;
 567     for (i = 0; i < _numSeqs; i++)
 568     {
 569         seqWeight[i] = temp;
 570     }
 571
 572     userParameters->setDistanceTree(false);
 573
 574     /* do the initial alignment.........  */
 575
 576     group.resize(_numSeqs + 1);
 577
 578     for (i = 1; i <= alnPtr->getProfile1NumSeqs(); ++i)
 579     {
 580         group[i] = 1;
 581     }
 582
 583     for (i = alnPtr->getProfile1NumSeqs() + 1; i <= _numSeqs; ++i)
 584     {
 585         group[i] = 2;
 586     }
 587
 588     entries = _numSeqs;
 589
 590     aligned = new int[_numSeqs + 1];
 591
 592     for (i = 1; i <= _numSeqs; i++)
 593     {
 594         aligned[i] = 1;
 595     }
 596
 597     alnPtr->addSeqWeight(&seqWeight);
 598
 599     score = alignAlg->profileAlign(alnPtr, distMat, &group, aligned);
 600     utilityObject->info("Sequences:%d      Score:%d", entries, score);
 601     delete [] aligned;
 602
 603     for (i = 1; i <= _numSeqs; i++)
 604     {
 605         for (j = i + 1; j <= _numSeqs; j++)
 606         {
 607             dscore = alnPtr->countid(i, j);
 608             (*distMat)(i, j) = ((double)100.0 - (double)dscore) / (double)100.0;
 609             (*distMat)(j, i) = (*distMat)(i, j);
 610         }
 611     }
 612     delete alignAlg;
 613     return (_numSeqs);
 614
 615 }
 616
 617
 618 /**
 619  *
 620  * @param alnPtr
 621  * @param distMat
 622  * @param p1TreeName
 623  * @param p2TreeName
 624  * @return
 625  */
 626 int MSA::doProfileAlign(Alignment* alnPtr, DistMatrix* distMat, vector<int>* prof1Weight, vector<int>* prof2Weight)
 627 {
 628     //Tree groupTree1, groupTree2;
 629     int i, j, sum, entries;
 630     int score;
 631     int *aligned;
 632     vector<int> group;
 633     int *maxid;
 634     //vector<int> prof1Weight, prof2Weight;
 635     int _profile1NumSeqs = alnPtr->getProfile1NumSeqs();
 636     int _numSeqs = alnPtr->getNumSeqs();
 637     vector<int> _seqWeight, _outputIndex;
 638
 639     utilityObject->info("Start of Multiple Alignment\n");
 640
 641     _seqWeight.resize(_numSeqs + 1);
 642     _outputIndex.resize(_numSeqs);
 643     ProfileAlignAlgorithm* alignAlgorithm = new MyersMillerProfileAlign;
 644
 645     // weight sequences with max percent identity with other profile
 646
 647     maxid = new int[_numSeqs + 1];
 648     for (i = 0; i < _profile1NumSeqs; i++)
 649     {
 650         maxid[i] = 0;
 651         for (j = _profile1NumSeqs + 1; j <= _numSeqs; j++)
 652         {
 653             if (maxid[i] < (*distMat)(i + 1, j))
 654             {
 655                 maxid[i] = static_cast<int>((*distMat)(i + 1, j)); // Mark change 17-5-07
 656             }
 657         }
 658         _seqWeight[i] = maxid[i] * (*prof1Weight)[i];
 659     }
 660
 661     for (i = _profile1NumSeqs; i < _numSeqs; i++)
 662     {
 663         maxid[i] =  - 1;
 664         for (j = 1; j <= _profile1NumSeqs; j++)
 665         {
 666             if (maxid[i] < (*distMat)(i + 1, j))
 667             {
 668                 maxid[i] = static_cast<int>((*distMat)(i + 1, j));// Mark change 17-5-07
 669             }
 670         }
 671         _seqWeight[i] = maxid[i] * (*prof2Weight)[i];
 672     }
 673     //
 674     // Normalise the weights, such that the sum of the weights = INT_SCALE_FACTOR
 675     //
 676
 677     sum = 0;
 678     for (j = 0; j < _numSeqs; j++)
 679     {
 680         sum += _seqWeight[j];
 681     }
 682
 683     if (sum == 0)
 684     {
 685         for (j = 0; j < _numSeqs; j++)
 686         {
 687             _seqWeight[j] = 1;
 688         }
 689
 690         sum = j;
 691     }
 692     for (j = 0; j < _numSeqs; j++)
 693     {
 694         _seqWeight[j] = (_seqWeight[j] * INT_SCALE_FACTOR) / sum;
 695         if (_seqWeight[j] < 1)
 696         {
 697             _seqWeight[j] = 1;
 698         }
 699     }
 700
 701     // do the alignment.........  /
 702
 703     utilityObject->info("Aligning...");
 704     group.resize(_numSeqs + 1);
 705
 706     for (i = 1; i <= _profile1NumSeqs; ++i)
 707     {
 708         group[i] = 1;
 709     }
 710
 711     for (i = _profile1NumSeqs + 1; i <= _numSeqs; ++i)
 712     {
 713         group[i] = 2;
 714     }
 715
 716     entries = _numSeqs;
 717
 718     aligned = new int[_numSeqs + 1];
 719     for (i = 1; i <= _numSeqs; i++)
 720     {
 721         aligned[i] = 1;
 722     }
 723
 724     alnPtr->addSeqWeight(&_seqWeight);
 725
 726     score = alignAlgorithm->profileAlign(alnPtr, distMat, &group, aligned);
 727
 728     utilityObject->info("Sequences:%d      Score:%d", entries, score);
 729
 730     for (i = 1; i <= _numSeqs; i++)
 731     {
 732         _outputIndex[i - 1] = i;
 733     }
 734
 735     alnPtr->addOutputIndex(&_outputIndex);
 736
 737     delete alignAlgorithm;
 738     delete [] aligned;
 739     delete [] maxid;
 740     return (_numSeqs);
 741     return 1;
 742 }
 743
 744 }