binaries/src/clustalw/src/pairwise/FastPairwiseAlign.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 <math.h>
  10 #include "FastPairwiseAlign.h"
  11
  12 namespace clustalw
  13 {
  14
  15 FastPairwiseAlign::FastPairwiseAlign()
  16 {
  17     _maxAlnLength = 0;
  18 }
  19
  20 void FastPairwiseAlign::pairwiseAlign(Alignment *alignPtr, DistMatrix *distMat, int iStart,
  21                                       int iEnd, int jStart, int jEnd)
  22 {
  23     try
  24     {
  25         if(distMat->getSize() != alignPtr->getNumSeqs() + 1)
  26         {
  27             cerr << "The distance matrix is not the right size!\n"
  28                  << "Need to terminate program.\n";
  29             exit(1);
  30         }
  31         if((iStart < 0) || (iEnd < iStart) || (jStart < 0) || (jEnd < jStart))
  32         {
  33             cout << "The range for pairwise Alignment is incorrect.\n"
  34                  << "Need to terminate program.\n";
  35             exit(1);
  36         }
  37
  38         int i, j, dsr;
  39         double calcScore;
  40         bool _DNAFlag = userParameters->getDNAFlag();
  41         _maxAlnLength = alignPtr->getMaxAlnLength();
  42         int num = (2 * _maxAlnLength) + 1;
  43         accum.ResizeRect(5, num);
  44
  45         displ.resize(num);
  46         slopes.resize(num);
  47         diagIndex.resize(num);
  48
  49         zza.resize(_maxAlnLength + 1);
  50         zzb.resize(_maxAlnLength + 1);
  51         zzc.resize(_maxAlnLength + 1);
  52         zzd.resize(_maxAlnLength + 1);
  53
  54         if (_DNAFlag)
  55         {
  56             userParameters->setDNAParams();
  57         }
  58         else
  59         {
  60             userParameters->setProtParams();
  61         }
  62
  63         cout << "\n\n";
  64
  65         for (i = iStart + 1; i <= iEnd; ++i)
  66         {
  67             const vector<int>* _seqIPtr = alignPtr->getSequence(i);
  68             int _seqILength = alignPtr->getSeqLength(i);
  69             if (_DNAFlag)
  70             {
  71                 makeNPtrs(zza, zzc, _seqIPtr, _seqILength);
  72             }
  73             else
  74             {
  75                 makePPtrs(zza, zzc, _seqIPtr, _seqILength);
  76             }
  77             double _score;
  78             for (j = jStart + 2 > i+1 ? jStart + 2 : i+1; j <= jEnd; ++j)
  79             {
  80                 const vector<int>* _seqJPtr = alignPtr->getSequence(j);
  81                 int _seqJLength = alignPtr->getSeqLength(j);
  82                 if (_DNAFlag)
  83                 {
  84                     makeNPtrs(zzb, zzd, _seqJPtr, _seqJLength);
  85                 }
  86                 else
  87                 {
  88                     makePPtrs(zzb, zzd, _seqJPtr, _seqJLength);
  89                 }
  90                 pairAlign(_seqIPtr, _seqILength, _seqJLength);
  91                 if (!maxSoFar)
  92                 {
  93                     calcScore = 0.0;
  94                 }
  95                 else
  96                 {
  97                     calcScore = (double)accum[0][maxSoFar];
  98                     if (userParameters->getPercent())
  99                     {
 100                         dsr = (_seqILength < _seqJLength) ? _seqILength
 101                             : _seqJLength;
 102                         calcScore = (calcScore / (double)dsr) *100.0;
 103                     }
 104                 }
 105                 _score = (100.0 - calcScore) / 100.0;
 106                 distMat->SetAt(i, j, _score);
 107                 //distMat->SetAt(j, i, _score); /* distMat symmetric, FS, 2009-04-06 */
 108
 109
 110                 if(userParameters->getDisplayInfo())
 111                 {
 112                     if (calcScore > 0.1)
 113                     {
 114                         utilityObject->info("Sequences (%d:%d) Aligned. Score: %lg",
 115                                             i, j, calcScore);
 116                     }
 117                     else
 118                     {
 119                         utilityObject->info("Sequences (%d:%d) Not Aligned", i, j);
 120                     }
 121                 }
 122             }
 123         }
 124         accum.clearArray();
 125         displ.clear();
 126         slopes.clear();
 127         diagIndex.clear();
 128
 129         zza.clear();
 130         zzb.clear();
 131         zzc.clear();
 132         zzd.clear();
 133     }
 134     catch(const exception& e)
 135     {
 136         cerr << "An exception has occured in the FastPairwiseAlign class.\n"
 137              << e.what() << "\n";
 138         exit(1);
 139     }
 140 }
 141
 142
 143 /*
 144  * Note: There is a problem with the treatment of DNA/RNA.
 145  * During file reading all residues are encoded as AminoAcids,
 146  * even before it has been established if they are AA or not.
 147  * This is bad and will have to be changed (later).
 148  * 'A' is assigned code 0, C is assigned 2, G = 6, T=18, U=19.
 149  * However, the fast alignment routines require that
 150  * A=0, C=1, G=2, T=U=3. In the best case the results of the
 151  * fast alignment (of DNA) will simply be meaningless, the
 152  * worst case is a core dump.
 153  * As a quick fix I implemented the following mask, that
 154  * (for DNA/RNA) translates 0->0, 2->1, 6->2, 18->3, 19->3.
 155  * This is awfull (it is not OO compliant) but it was quick,
 156  * uses (much) less memory than a second DNA array, and is
 157  * (much) faster than calling a translation function.
 158  * Ideally this will be removed, but this requires changes
 159  * to (i) the sequence encoding during file-reading AND
 160  * (ii) all the DNA substitution matrices. (Fabian, 2009-02-25)
 161  *
 162  *                A  B  C  D  E  F  G  H  I  K  L  M  0  N  P  Q  R  S
 163  *                T  U  W  X  Y  Z (goodmeasuregoodmeasuregood) */
 164 int ziAA2DNA[] = {0,-1, 1,-1,-1,-1, 2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
 165                   3, 3,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
 166 const int ciMaxResID = 3;
 167
 168
 169
 170
 171 void FastPairwiseAlign::pairAlign(const vector<int>* seq, int l1, int l2)
 172 {
 173     int pot[8], i, j, l, m, limit, pos, tl1, vn1, vn2, flen, osptr, fs;
 174     int tv1, tv2, encrypt, subt1, subt2, rmndr;
 175     bool flag;
 176     int residue;
 177     bool _DNAFlag = userParameters->getDNAFlag();
 178     int _ktup = userParameters->getKtup();
 179     int _maxAA = userParameters->getMaxAA();
 180     int _windowGap = userParameters->getWindowGap();
 181     int _window = userParameters->getWindow();
 182     int _signif = userParameters->getSignif();
 183
 184     if (_DNAFlag)
 185     {
 186         for (i = 1; i <= _ktup; ++i)
 187         {
 188             pot[i] = (int)pow((double)4, (double)(i - 1));
 189         }
 190         limit = (int)pow((double)4, (double)_ktup);
 191     }
 192     else
 193     {
 194         for (i = 1; i <= _ktup; i++)
 195         {
 196             pot[i] = (int)pow((double)(_maxAA + 1), (double)(i - 1));
 197         }
 198         limit = (int)pow((double)(_maxAA + 1), (double)_ktup);
 199     }
 200     tl1 = (l1 + l2) - 1;
 201
 202     for (i = 1; i <= tl1; ++i)
 203     {
 204         slopes[i] = displ[i] = 0;
 205         diagIndex[i] = i;
 206     }
 207
 208
 209     // increment diagonal score for each k_tuple match
 210
 211     for (i = 1; i <= limit; ++i)
 212     {
 213         vn1 = zzc[i];
 214         while (true)
 215         {
 216             if (!vn1)
 217             {
 218                 break;
 219             }
 220             vn2 = zzd[i];
 221             while (vn2 != 0)
 222             {
 223                 osptr = vn1 - vn2 + l2;
 224                 ++displ[osptr];
 225                 vn2 = zzb[vn2];
 226             }
 227             vn1 = zza[vn1];
 228         }
 229     }
 230
 231     // choose the top SIGNIF diagonals
 232
 233     desQuickSort(displ, diagIndex, tl1);
 234
 235     j = tl1 - _signif + 1;
 236     if (j < 1)
 237     {
 238         j = 1;
 239     }
 240
 241     // flag all diagonals within WINDOW of a top diagonal
 242
 243     for (i = tl1; i >= j; i--)
 244     if (displ[i] > 0)
 245     {
 246         pos = diagIndex[i];
 247         l = (1 > pos - _window) ? 1 : pos - _window;
 248         m = (tl1 < pos + _window) ? tl1 : pos + _window;
 249         for (; l <= m; l++)
 250         {
 251             slopes[l] = 1;
 252         }
 253     }
 254
 255     for (i = 1; i <= tl1; i++)
 256     {
 257         displ[i] = 0;
 258     }
 259
 260
 261     currFrag = maxSoFar = 0;
 262
 263     for (i = 1; i <= (l1 - _ktup + 1); ++i)
 264     {
 265         encrypt = flag = 0;
 266         if (_DNAFlag){
 267           for (j = 1; j <= _ktup; ++j)
 268             {
 269               residue = ziAA2DNA[(*seq)[i + j - 1]];
 270               if ((residue < 0) || (residue > ciMaxResID))
 271                 {
 272                   flag = true;
 273                   break;
 274                 }
 275               encrypt += ((residue) * pot[j]);
 276             }
 277         }
 278         else {
 279           for (j = 1; j <= _ktup; ++j)
 280             {
 281               residue = (*seq)[i + j - 1];
 282               if ((residue < 0) || (residue > _maxAA))
 283                 {
 284                   flag = true;
 285                   break;
 286                 }
 287               encrypt += ((residue) * pot[j]);
 288             }
 289         }
 290
 291         if (flag)
 292         {
 293             continue;
 294         }
 295         ++encrypt;
 296
 297         vn2 = zzd[encrypt];
 298
 299         flag = false;
 300         while (true)
 301         {
 302             if (!vn2)
 303             {
 304                 flag = true;
 305                 break;
 306             }
 307             osptr = i - vn2 + l2;
 308             if (slopes[osptr] != 1)
 309             {
 310                 vn2 = zzb[vn2];
 311                 continue;
 312             }
 313             flen = 0;
 314             fs = _ktup;
 315             next = maxSoFar;
 316
 317
 318             //  A-loop
 319
 320
 321             while (true)
 322             {
 323                 if (!next)
 324                 {
 325                     ++currFrag;
 326                     if (currFrag >= 2 * _maxAlnLength)
 327                     {
 328                         utilityObject->info("(Partial alignment)");
 329                         vatend = 1;
 330                         return ;
 331                     }
 332                     displ[osptr] = currFrag;
 333                     putFrag(fs, i, vn2, flen);
 334                 }
 335                 else
 336                 {
 337                     tv1 = accum[1][next];
 338                     tv2 = accum[2][next];
 339                     if (fragRelPos(i, vn2, tv1, tv2))
 340                     {
 341                         if (i - vn2 == accum[1][next] - accum[2][next])
 342                         {
 343                             if (i > accum[1][next] + (_ktup - 1))
 344                             {
 345                                 fs = accum[0][next] + _ktup;
 346                             }
 347                             else
 348                             {
 349                                 rmndr = i - accum[1][next];
 350                                 fs = accum[0][next] + rmndr;
 351                             }
 352                             flen = next;
 353                             next = 0;
 354                             continue;
 355                         }
 356                         else
 357                         {
 358                             if (displ[osptr] == 0)
 359                             {
 360                                 subt1 = _ktup;
 361                             }
 362                             else
 363                             {
 364                                 if (i > accum[1][displ[osptr]] + (_ktup - 1))
 365                                 {
 366                                     subt1 = accum[0][displ[osptr]] + _ktup;
 367                                 }
 368                                 else
 369                                 {
 370                                     rmndr = i - accum[1][displ[osptr]];
 371                                     subt1 = accum[0][displ[osptr]] + rmndr;
 372                                 }
 373                             }
 374                             subt2 = accum[0][next] - _windowGap + _ktup;
 375                             if (subt2 > subt1)
 376                             {
 377                                 flen = next;
 378                                 fs = subt2;
 379                             }
 380                             else
 381                             {
 382                                 flen = displ[osptr];
 383                                 fs = subt1;
 384                             }
 385                             next = 0;
 386                             continue;
 387                         }
 388                     }
 389                     else
 390                     {
 391                         next = accum[4][next];
 392                         continue;
 393                     }
 394                 }
 395                 break;
 396             }
 397             // End of Aloop
 398             vn2 = zzb[vn2];
 399         }
 400     }
 401     vatend = 0;
 402 }
 403
 404 void FastPairwiseAlign::makePPtrs(vector<int>& tptr, vector<int>& pl, const vector<int>* seq, int length)
 405 {
 406     int a[10];
 407     int i, j, limit, code;
 408     bool flag;
 409     int residue;
 410     int _ktup = userParameters->getKtup();
 411     int _maxAA = userParameters->getMaxAA();
 412
 413     for (i = 1; i <= _ktup; i++)
 414     {
 415         a[i] = (int)pow((double)(_maxAA + 1), (double)(i - 1));
 416     }
 417
 418     limit = (int)pow((double)(_maxAA + 1), (double)_ktup);
 419     if(limit >= (int)pl.size())
 420     {
 421         pl.resize(limit + 1);
 422     }
 423     if(length >= (int)tptr.size())
 424     {
 425         tptr.resize(length + 1);
 426     }
 427
 428     for (i = 1; i <= limit; ++i)
 429     {
 430         pl[i] = 0;
 431     }
 432     for (i = 1; i <= length; ++i) // NOTE changed this
 433     {
 434         tptr[i] = 0;
 435     }
 436
 437     for (i = 1; i <= (length - _ktup + 1); ++i)
 438     {
 439         code = 0;
 440         flag = false;
 441         for (j = 1; j <= _ktup; ++j)
 442         {
 443             residue = (*seq)[i + j - 1];
 444             if ((residue < 0) || (residue > _maxAA))
 445             {
 446                 flag = true;
 447                 break;
 448             }
 449             code += ((residue) * a[j]);
 450         }
 451         if (flag)
 452         {
 453             continue;
 454         }
 455         ++code;
 456         if (pl[code] != 0)
 457         {
 458             tptr[i] = pl[code];
 459         }
 460         pl[code] = i; //Ktuple code is at position i in sequence
 461     }
 462 }
 463
 464 void FastPairwiseAlign::makeNPtrs(vector<int>& tptr,vector<int>& pl, const vector<int>* seq, int length)
 465 {
 466     int pot[] =
 467     {
 468         0, 1, 4, 16, 64, 256, 1024, 4096
 469     };
 470     int i, j, limit, code;
 471     bool flag;
 472     int residue;
 473     int _ktup = userParameters->getKtup();
 474
 475     limit = (int)pow((double)4, (double)_ktup);
 476
 477     if(limit >= (int)pl.size())
 478     {
 479         pl.resize(limit + 1);
 480     }
 481     if(length >= (int)tptr.size())
 482     {
 483         tptr.resize(length + 1);
 484     }
 485
 486     for (i = 1; i <= limit; ++i)
 487     {
 488         pl[i] = 0;
 489     }
 490     for (i = 1; i <= length; ++i)
 491     {
 492         tptr[i] = 0;
 493     }
 494
 495     for (i = 1; i <= length - _ktup + 1; ++i)
 496     {
 497         code = 0;
 498         flag = false;
 499         for (j = 1; j <= _ktup; ++j)
 500           {
 501             residue = ziAA2DNA[(*seq)[i + j - 1]];
 502             if ((residue < 0) || (residue > ciMaxResID))
 503               {
 504                 flag = true;
 505                 break;
 506               }
 507             code += ((residue) * pot[j]);
 508           }
 509         if (flag)
 510         {
 511             continue;
 512         }
 513         ++code;
 514         if (pl[code] != 0)
 515         {
 516             tptr[i] = pl[code];
 517         }
 518         pl[code] = i;
 519     }
 520
 521 }
 522
 523 void FastPairwiseAlign::putFrag(int fs, int v1, int v2, int flen)
 524 {
 525     int end;
 526     accum[0][currFrag] = fs;
 527     accum[1][currFrag] = v1;
 528     accum[2][currFrag] = v2;
 529     accum[3][currFrag] = flen;
 530
 531     if (!maxSoFar)
 532     {
 533         maxSoFar = 1;
 534         accum[4][currFrag] = 0;
 535         return ;
 536     }
 537
 538     if (fs >= accum[0][maxSoFar])
 539     {
 540         accum[4][currFrag] = maxSoFar;
 541         maxSoFar = currFrag;
 542         return ;
 543     }
 544     else
 545     {
 546         next = maxSoFar;
 547         while (true)
 548         {
 549             end = next;
 550             next = accum[4][next];
 551             if (fs >= accum[0][next])
 552             {
 553                 break;
 554             }
 555         }
 556         accum[4][currFrag] = next;
 557         accum[4][end] = currFrag;
 558     }
 559 }
 560
 561 inline int FastPairwiseAlign::fragRelPos(int a1, int b1, int a2, int b2)
 562 {
 563     int ret;
 564     int _ktup = userParameters->getKtup();
 565
 566     ret = false;
 567     if (a1 - b1 == a2 - b2)
 568     {
 569         if (a2 < a1)
 570         {
 571             ret = true;
 572         }
 573     }
 574     else
 575     {
 576         if (a2 + _ktup - 1 < a1 && b2 + _ktup - 1 < b1)
 577         {
 578             ret = true;
 579         }
 580     }
 581     return ret;
 582 }
 583
 584 void FastPairwiseAlign::desQuickSort(vector<int>& array1, vector<int>& array2, int arraySize)
 585 {
 586     // Quicksort routine, adapted from chapter 4, page 115 of software tools
 587     // by Kernighan and Plauger, (1986)
 588     // Sort the elements of array1 and sort the
 589     // elements of array2 accordingly
 590     int temp1, temp2;
 591     int p, pivlin;
 592     int i, j;
 593     int lst[50], ust[50]; // the maximum no. of elements must be
 594                             // < log(base2) of 50
 595
 596     lst[1] = 1;
 597     ust[1] = arraySize - 1;
 598     p = 1;
 599
 600     while (p > 0)
 601     {
 602         if (lst[p] >= ust[p])
 603         {
 604             p--;
 605         }
 606         else
 607         {
 608             i = lst[p] - 1;
 609             j = ust[p];
 610             pivlin = array1[j];
 611             while (i < j)
 612             {
 613                 for (i = i + 1; array1[i] < pivlin; i++)
 614                     ;
 615                 for (j = j - 1; j > i; j--)
 616                     if (array1[j] <= pivlin)
 617                     {
 618                         break;
 619                     }
 620                 if (i < j)
 621                 {
 622                     temp1 = array1[i];
 623                     array1[i] = array1[j];
 624                     array1[j] = temp1;
 625
 626                     temp2 = array2[i];
 627                     array2[i] = array2[j];
 628                     array2[j] = temp2;
 629                 }
 630             }
 631
 632             j = ust[p];
 633
 634             temp1 = array1[i];
 635             array1[i] = array1[j];
 636             array1[j] = temp1;
 637
 638             temp2 = array2[i];
 639             array2[i] = array2[j];
 640             array2[j] = temp2;
 641
 642             if (i - lst[p] < ust[p] - i)
 643             {
 644                 lst[p + 1] = lst[p];
 645                 ust[p + 1] = i - 1;
 646                 lst[p] = i + 1;
 647             }
 648             else
 649             {
 650                 lst[p + 1] = i + 1;
 651                 ust[p + 1] = ust[p];
 652                 ust[p] = i - 1;
 653             }
 654             p = p + 1;
 655         }
 656     }
 657     return ;
 658
 659 }
 660
 661 }