binaries/src/clustalw/src/alignment/Alignment.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 /**
   7  * Changes:
   8  *
   9  * 16-02-07,Nigel Brown(EMBL): Added friend NameIterator to allow a caller to
  10  * process the name vector.
  11  *
  12  * 05-03-07,Nigel Brown(EMBL): modified searchForString() to skip over gaps in
  13  * sequence while matching.
  14  *
  15  * 26-03-07,Nigel Brown(EMBL): suppressed error message box for name clashes;
  16  * added testUniqueNames() predicate, which compares new sequence names with
  17  * those in the alignment vector BEFORE appending them.
  18  */
  19 #ifdef HAVE_CONFIG_H
  20     #include "config.h"
  21 #endif
  22 #include <exception>
  23 #include <cmath>
  24 #include <sstream>
  25 #include "Alignment.h"
  26
  27 using namespace std;
  28
  29 namespace clustalw
  30 {
  31
  32 Alignment::Alignment()
  33  : gapPos1(userParameters->getGapPos1()),
  34    gapPos2(userParameters->getGapPos2())
  35 {
  36     maxNames = 0;
  37     numSeqs = 0;
  38     maxAlignmentLength = 0;
  39     lengthLongestSequence = 0;
  40     profile1NumSeqs = 0;
  41 }
  42
  43 // Andreas Wilm (UCD): shouldn't resetProfile1 and resetProfile2 be merged?
  44 /** remove gaps from older alignments (code = gap_pos1) */
  45 void Alignment::resetProfile1()
  46 {
  47     register int sl;                /* which have  code = gap_pos2  */
  48     int i,j;
  49     int _gapPos1 = userParameters->getGapPos1();
  50     int _gapPos2 = userParameters->getGapPos2();
  51     bool _resetAlignNew = userParameters->getResetAlignmentsNew();
  52     bool _resetAlignAll = userParameters->getResetAlignmentsAll();
  53
  54     if (userParameters->getStructPenalties1() != NONE)
  55     {
  56         sl = 0;
  57         for (j = 0; j < (int)gapPenaltyMask1.size(); ++j)
  58         {
  59             if (gapPenaltyMask1[j] == _gapPos1 && (_resetAlignNew ||
  60                            _resetAlignAll))
  61             {
  62                 continue;
  63             }
  64             if (gapPenaltyMask1[j] == _gapPos2 && (_resetAlignAll))
  65             {
  66                 continue;
  67             }
  68             gapPenaltyMask1[sl] = gapPenaltyMask1[j];
  69             ++sl;
  70         }
  71     }
  72
  73     if (userParameters->getStructPenalties1() == SECST)
  74     {
  75         sl = 0;
  76         for (j = 0; j < (int)secStructMask1.size(); ++j)
  77         {
  78             if (secStructMask1[j] == _gapPos1 && (_resetAlignNew ||
  79                           _resetAlignAll))
  80             {
  81                 continue;
  82             }
  83             if (secStructMask1[j] == _gapPos2 && (_resetAlignAll))
  84             {
  85                 continue;
  86             }
  87             secStructMask1[sl] = secStructMask1[j];
  88             ++sl;
  89         }
  90     }
  91
  92     for(i = 1; i <= profile1NumSeqs; ++i)
  93     {
  94         sl = 0;
  95         for(j = 1; j <= getSeqLength(i); ++j)
  96         {
  97             if(seqArray[i][j] == _gapPos1 && (_resetAlignNew || _resetAlignAll))
  98             {
  99                 continue;
 100             }
 101             if(seqArray[i][j] == _gapPos2 && (_resetAlignAll))
 102             {
 103                 continue;
 104             }
 105             ++sl;
 106             seqArray[i][sl] = seqArray[i][j];
 107         }
 108
 109         // Andreas Wilm (UCD): added 2008-03-07:
 110         // Remove excess bit at end of sequence
 111         int numExtraElements = seqArray[i].size() - 1 - sl;
 112         for(int k = 0; k < numExtraElements; k++)
 113         {
 114             seqArray[i].pop_back();
 115         }
 116     }
 117 }
 118
 119 // Andreas Wilm (UCD): shouldn't resetProfile1 and resetProfile2 be merged?
 120 void Alignment::resetProfile2()
 121 {
 122     register int sl;                /* which have  code = gap_pos2  */
 123     int i, j;
 124     int _gapPos1 = userParameters->getGapPos1();
 125     int _gapPos2 = userParameters->getGapPos2();
 126     bool _resetAlignNew = userParameters->getResetAlignmentsNew();
 127     bool _resetAlignAll = userParameters->getResetAlignmentsAll();
 128     int _profile1NumSeqs = profile1NumSeqs;
 129
 130
 131     if (userParameters->getStructPenalties2() != NONE)
 132     {
 133         sl = 0;
 134         for (j = 0; j < (int)gapPenaltyMask2.size(); ++j)
 135         {
 136             if (gapPenaltyMask2[j] == _gapPos1 && (_resetAlignNew || _resetAlignAll))
 137             {
 138                 continue;
 139             }
 140             if (gapPenaltyMask2[j] == _gapPos2 && (_resetAlignAll))
 141             {
 142                 continue;
 143             }
 144             gapPenaltyMask2[sl] = gapPenaltyMask2[j];
 145             ++sl;
 146         }
 147     }
 148
 149     if (userParameters->getStructPenalties2() == SECST)
 150     {
 151         sl = 0;
 152         for (j = 0; j < (int)secStructMask2.size(); ++j)
 153         {
 154             if (secStructMask2[j] == _gapPos1 && (_resetAlignNew ||
 155                          _resetAlignAll))
 156             {
 157                 continue;
 158             }
 159             if (secStructMask2[j] == _gapPos2 && (_resetAlignAll))
 160             {
 161                 continue;
 162             }
 163             secStructMask2[sl] = secStructMask2[j];
 164             ++sl;
 165         }
 166     }
 167
 168     for(i = _profile1NumSeqs + 1; i <= numSeqs; ++i)
 169     {
 170         sl = 0;
 171         for(j = 1; j <= getSeqLength(i); ++j)
 172         {
 173             if(seqArray[i][j] == _gapPos1 && (_resetAlignNew || _resetAlignAll))
 174             {
 175                 continue;
 176             }
 177             if(seqArray[i][j] == _gapPos2 && (_resetAlignAll))
 178             {
 179                 continue;
 180             }
 181             ++sl;
 182             seqArray[i][sl] = seqArray[i][j];
 183         }
 184
 185
 186         // Andreas Wilm (UCD) added 2008-03-07:
 187         // Remove excess bit at end of sequence
 188         int numExtraElements = seqArray[i].size() - 1 - sl;
 189         for(int k = 0; k < numExtraElements; k++)
 190         {
 191             seqArray[i].pop_back();
 192         }
 193     }
 194
 195 }
 196
 197 void Alignment::resetAllSeqWeights()
 198 {
 199     seqWeight.clear();
 200     seqWeight.resize(numSeqs + 1, 100);
 201 }
 202
 203 /*
 204  * The function addOutputIndex must check that the number of outputindex
 205  * Not sure what to do with this one.
 206  */
 207 bool Alignment::addOutputIndex(vector<int>* outputIndexToAdd)
 208 {
 209     // First need to clear the old outputIndex
 210     outputIndex.clear();
 211     // Check if the size is correct
 212     if((int)outputIndexToAdd->size() == numSeqs)
 213     {
 214         // Add the output index
 215         outputIndex = *outputIndexToAdd;
 216         return true;
 217     }
 218     else
 219     {
 220         clearAlignment();
 221         return false; // Could not add them
 222     }
 223
 224 }
 225
 226 /*
 227  * The function appendOutputIndex is used when we are appending the outputIndex
 228  * to the current one. We do this when we are adding a profile2.
 229  */
 230 bool Alignment::appendOutputIndex(vector<int>* outputIndexToAppend)
 231 {
 232     // Check if the size is correct
 233     if((int)(outputIndexToAppend->size() + outputIndex.size()) == numSeqs)
 234     {
 235         // Append the outputIndex
 236         vector<int>::iterator outIndexIter = outputIndexToAppend->begin();
 237         while(outIndexIter != outputIndexToAppend->end())
 238         {
 239             outputIndex.push_back(*outIndexIter);
 240             outIndexIter++;
 241         }
 242         if((int)outputIndex.size() == numSeqs)
 243         {
 244             return true;
 245         }
 246         else
 247         {
 248             clearAlignment();
 249             cerr << "There is a problem with adding the sequences\n";
 250             return false;
 251         }
 252     }
 253     else
 254     {
 255         clearAlignment();
 256         return false;
 257     }
 258 }
 259
 260 /*
 261  * The function addSequences takes a vector of sequences and adds it to the Alignment.
 262  * This is used if there is nothing in memory that we need.
 263  *
 264  */
 265 void Alignment::addSequences(vector<Sequence>* seqVector)
 266 {
 267     clearAlignment();
 268
 269     numSeqs = seqVector->size();
 270     vector<int> emptyVec;
 271
 272     // NOTE push dummy sequences on first!
 273     /********************************************************
 274      * It was decided to stay with the seqs and seqArray    *
 275      * index begining at 1. This is because it is difficult *
 276      * to change in the code, so we push dummy seqs on      *
 277      ********************************************************/
 278     seqArray.push_back(emptyVec); // EMPTY sequence
 279     names.push_back(string(""));
 280     titles.push_back(string(""));
 281     sequenceIds.push_back(0);
 282
 283     addSequencesToVector(seqVector);
 284
 285     calculateMaxLengths();
 286     seqWeight.resize(numSeqs + 1, 100);
 287 }
 288
 289 /*
 290  * The function appendSequences is used when we have a profile2.
 291  *
 292  */
 293 void Alignment::appendSequences(vector<Sequence>* seqVector)
 294 {
 295     numSeqs += seqVector->size(); // Add to the number we already have!
 296     addSequencesToVector(seqVector);
 297     seqWeight.clear();
 298     seqWeight.resize(numSeqs + 1, 100);
 299     calculateMaxLengths();
 300 }
 301
 302 void Alignment::pasteSequencesIntoPosition(vector<Sequence>* seqVector, int pos,
 303                                            bool explicitPasteToProfile2)
 304 {
 305     SeqArray::iterator seqArrayIterator;
 306     vector<string>::iterator namesIterator;
 307     vector<string>::iterator titlesIterator;
 308     vector<unsigned long>::iterator sequenceIdsIterator;
 309
 310     int profNum = userParameters->getProfileNum();
 311     int numSeqsToInsert = seqVector->size();
 312     if(numSeqsToInsert == 0 || pos < 0)
 313     {
 314         return;
 315     }
 316     if(pos == numSeqs)
 317     {
 318         seqArrayIterator = seqArray.end();
 319         namesIterator = names.end();
 320         titlesIterator = titles.end();
 321         sequenceIdsIterator = sequenceIds.end();
 322     }
 323     else
 324     {
 325         seqArrayIterator = seqArray.begin() + pos + 1;
 326         namesIterator = names.begin() + pos + 1;
 327         titlesIterator = titles.begin() + pos + 1;
 328         sequenceIdsIterator = sequenceIds.begin() + pos + 1;
 329     }
 330
 331     int prof1NumSeqs = profile1NumSeqs;
 332
 333     for(int i = numSeqsToInsert - 1; i >= 0; i--)
 334     {
 335         seqArray.insert(seqArrayIterator, *(*seqVector)[i].getSequence());
 336         names.insert(namesIterator, (*seqVector)[i].getName());
 337         titles.insert(titlesIterator, (*seqVector)[i].getTitle());
 338         sequenceIds.insert(sequenceIdsIterator, (*seqVector)[i].getIdentifier());
 339
 340         numSeqs++;
 341         if(profNum != 0 && !explicitPasteToProfile2 && pos <= prof1NumSeqs)
 342         {
 343             prof1NumSeqs++;
 344         }
 345     }
 346
 347     if(profNum != 0 && pos <= prof1NumSeqs)
 348     {
 349         profile1NumSeqs = prof1NumSeqs;
 350     }
 351
 352     resetAllSeqWeights();
 353     setDefaultOutputIndex();
 354 }
 355
 356 void Alignment::debugPrintAllNames()
 357 {
 358     vector<string>::iterator nameIter = names.begin();
 359     while(nameIter != names.end())
 360     {
 361         cout << *nameIter << endl;
 362         nameIter++;
 363     }
 364 }
 365
 366 void Alignment::NameIterator::begin(Alignment *a) {
 367     //cout << "begin()\n";
 368     if (a) {
 369         alignment = a;
 370         i = alignment->names.begin();
 371     }
 372 }
 373
 374 const string Alignment::NameIterator::next() {
 375     //cout << "next()\n";
 376     if (!alignment)
 377         return string();
 378     if (i == alignment->names.end())
 379         return string();
 380     return *i++;
 381 }
 382
 383 bool Alignment::NameIterator::end() {
 384     //cout << "end()\n";
 385     if (!alignment)
 386         return true;
 387     if (i == alignment->names.end())
 388         return true;
 389     return false;
 390 }
 391
 392 // 26-03-03,nige: test before appending seqVector
 393 bool Alignment::testUniqueNames(vector<Sequence>* seqVector, string *offendingSeq)
 394 {
 395     vector<string>::iterator   oldName;
 396     vector<Sequence>::iterator newName;
 397     bool unique = true;
 398
 399     //iterate over new candidate names
 400     for (newName = seqVector->begin(); unique && newName != seqVector->end(); newName++) {
 401         //iterate over old stored names
 402         for (oldName = names.begin() + 1; unique && oldName != names.end(); oldName++) {
 403             if (*oldName == newName->getName()) {
 404                 offendingSeq->assign(*oldName);
 405                 unique = false;
 406             }
 407         }
 408     }
 409     return unique;
 410 }
 411
 412
 413 /*
 414  * The function addSequencesToVector is used to add sequences to our seqVector
 415  * It is used by both addSequences and appendSequences. This is to reduce code
 416  * duplication.
 417  */
 418 void Alignment::addSequencesToVector(vector<Sequence>* seqVector)
 419 {
 420     std::vector<Sequence>::iterator itSeq;
 421
 422     for(itSeq = seqVector->begin(); itSeq != seqVector->end(); ++itSeq)
 423     {
 424         seqArray.push_back(*(*itSeq).getSequence());
 425         names.push_back((*itSeq).getName());
 426         titles.push_back((*itSeq).getTitle());
 427         sequenceIds.push_back((*itSeq).getIdentifier());
 428     }
 429
 430     if(!(((int)seqArray.size() == numSeqs + 1) && ((int)names.size() == numSeqs + 1)
 431           && ((int)titles.size() == numSeqs + 1) && ((int)sequenceIds.size() == numSeqs + 1)))
 432     {
 433         cerr << "There has been an error adding the sequences to Alignment.\n"
 434              << "Must terminate the program. EaddSequencesrror occured in addSequences.\n";
 435         exit(1);
 436     }
 437 }
 438
 439 void Alignment::clearAlignment()
 440 {
 441     // Erase all the elements from the vector!
 442     clearSeqArray();
 443     names.clear();
 444     titles.clear();
 445     outputIndex.clear();
 446     sequenceIds.clear();
 447     clearSecStruct1();
 448     clearSecStruct2();
 449     seqWeight.clear();
 450     maxNames = 0;
 451     numSeqs = 0;
 452     maxAlignmentLength = 0;
 453     lengthLongestSequence = 0;
 454     userParameters->setProfileNum(0);
 455     userParameters->setProfile1Empty(true);
 456     userParameters->setProfile2Empty(true);
 457 }
 458
 459 void Alignment::clearSeqArray()
 460 {
 461     for(int i = 0; i < (int)seqArray.size(); i++)
 462     {
 463         seqArray[i].clear();
 464     }
 465     seqArray.clear();
 466 }
 467
 468 void Alignment::clearSecStruct1()
 469 {
 470     gapPenaltyMask1.clear();
 471     secStructMask1.clear();
 472     secStructName1 = "";
 473 }
 474
 475 void Alignment::clearSecStruct2()
 476 {
 477     gapPenaltyMask2.clear();
 478     secStructMask2.clear();
 479     secStructName2 = "";
 480 }
 481
 482 /*
 483  * The function alignScore is used to score the alignment!
 484  * This is used by other classes to see what score the alignment has gotten.
 485  */
 486 int Alignment::alignScore(void)
 487 {
 488     int maxRes;
 489     int seq1, seq2, res1, res2;
 490     int ngaps;
 491     int i, len1, len2;
 492     int score;
 493     int _maxAA = userParameters->getMaxAA();
 494     float _gapOpen = userParameters->getGapOpen();
 495     int matrix[NUMRES][NUMRES];
 496
 497     //
 498     // calculate an overall score for the alignment by summing the
 499     // scores for each pairwise alignment
 500     //
 501     maxRes = subMatrix->getAlnScoreMatrix(matrix);
 502     if (maxRes == 0)
 503     {
 504         utilityObject->error("Matrix for alignment scoring not found\n");
 505         return 0;
 506     }
 507
 508     score = 0;
 509     for (seq1 = 1; seq1 <= numSeqs; seq1++)
 510     {
 511         for (seq2 = 1; seq2 < seq1; seq2++)
 512         {
 513             len1 = seqArray[seq1].size() - 1;
 514             len2 = seqArray[seq2].size() - 1;
 515             for (i = 1; i < len1 && i < len2; i++)
 516             {
 517                 res1 = seqArray[seq1][i];
 518                 res2 = seqArray[seq2][i];
 519
 520                 if ((res1 >= 0) && (res1 <= _maxAA) && (res2 >= 0) && (res2 <= _maxAA))
 521                 {
 522                     score += matrix[res1][res2];
 523                 }
 524             }
 525
 526             ngaps = countGaps(seq1, seq2, len1);
 527
 528             score = static_cast<int>(score - (100 * _gapOpen * ngaps)); // Mark change 17-5-07
 529
 530         }
 531     }
 532
 533     score /= 100;
 534
 535     utilityObject->info("Alignment Score %d\n", score);
 536     return score;
 537 }
 538
 539 int Alignment::countGaps(int seq1, int seq2, int len)
 540 {
 541     int i, g;
 542     int q, r;//,  *Q,  *R;
 543
 544     vector<int> Q, R;
 545     Q.resize(len + 2, 0);
 546     R.resize(len + 2, 0);
 547
 548     int _maxAA = userParameters->getMaxAA();
 549
 550     try
 551     {
 552         Q[0] = R[0] = g = 0;
 553
 554         for (i = 1; i < len; i++)
 555         {
 556             if (seqArray[seq1][i] > _maxAA)
 557             {
 558                 q = 1;
 559             }
 560             else
 561             {
 562                 q = 0;
 563             }
 564
 565             if (seqArray[seq2][i] > _maxAA)
 566             {
 567                 r = 1;
 568             }
 569             else
 570             {
 571                 r = 0;
 572             }
 573
 574             // NOTE I havent a clue what this does!!!!
 575             if (((Q[i - 1] <= R[i - 1]) && (q != 0) && (1-r != 0)) ||
 576                 ((Q[i - 1] >= R[i - 1]) && (1-q != 0) && (r != 0)))
 577             {
 578                 g += 1;
 579             }
 580
 581             if (q != 0)
 582             {
 583                 Q[i] = Q[i - 1] + 1;
 584             }
 585             else
 586             {
 587                 Q[i] = 0;
 588             }
 589
 590             if (r != 0)
 591             {
 592                 R[i] = R[i - 1] + 1;
 593             }
 594             else
 595             {
 596                 R[i] = 0;
 597             }
 598         }
 599     }
 600     catch(const exception &ex)
 601     {
 602         cerr << ex.what() << endl;
 603         cerr << "Terminating program. Cannot continue. Function = countGaps\n";
 604         exit(1);
 605     }
 606
 607     return (g);
 608 }
 609
 610 void Alignment::resetAlign()
 611 {
 612     /* remove gaps from older alignments (code = gap_pos1) EXCEPT for
 613        gaps that were INPUT with the seqs.  which have code =
 614        gap_pos2
 615     */
 616
 617     register int sl;
 618     int i, j;
 619     int _gapPos1 = userParameters->getGapPos1();
 620     int _gapPos2 = userParameters->getGapPos2();
 621     bool _resetAlignNew = userParameters->getResetAlignmentsNew();
 622     bool _resetAlignAll = userParameters->getResetAlignmentsAll();
 623
 624
 625     for(i = 1; i <= numSeqs;++i)
 626     {
 627         sl = 0;
 628         for(j = 1; j <= getSeqLength(i); ++j)
 629         {
 630             if(seqArray[i][j] == _gapPos1 && ( _resetAlignNew || _resetAlignAll))
 631             {
 632                 continue;
 633             }
 634             if(seqArray[i][j] == _gapPos2 && (_resetAlignAll))
 635             {
 636                 continue;
 637             }
 638             ++sl;
 639             seqArray[i][sl] = seqArray[i][j];
 640         }
 641
 642
 643         // Andreas Wilm (UCD) added 2008-03-07:
 644         // Remove excess bit at end of sequence
 645         int numExtraElements = seqArray[i].size() - 1 - sl;
 646         for(int k = 0; k < numExtraElements; k++)
 647         {
 648             seqArray[i].pop_back();
 649         }
 650     }
 651 }
 652
 653 void Alignment::fixGaps()
 654 {
 655     int i,j;
 656
 657     if (userParameters->getStructPenalties1() != NONE)
 658     {
 659         for (j = 0; j < getSeqLength(1); ++j)
 660         {
 661             if (gapPenaltyMask1[j] == userParameters->getGapPos1())
 662             {
 663                 gapPenaltyMask1[j] = userParameters->getGapPos2();
 664             }
 665         }
 666     }
 667
 668     if (userParameters->getStructPenalties1() == SECST)
 669     {
 670         for (j = 0; j < getSeqLength(1); ++j)
 671         {
 672             if (secStructMask1[j] == userParameters->getGapPos1())
 673             {
 674                 secStructMask1[j] = userParameters->getGapPos2();
 675             }
 676         }
 677     }
 678
 679     for(i = 1; i <= numSeqs; ++i)
 680     {
 681         for(j = 1; j <= getSeqLength(i); ++j)
 682         {
 683             if(seqArray[i][j] == userParameters->getGapPos1())
 684             {
 685                 seqArray[i][j] = userParameters->getGapPos2();
 686             }
 687         }
 688     }
 689
 690 }
 691
 692 float Alignment::countid(int s1, int s2)
 693 {
 694     int c1, c2;
 695     int i;
 696     int count, total;
 697     float score;
 698     int shorterSeqLength = (getSeqLength(s1) < getSeqLength(s2)) ? getSeqLength(s1) : getSeqLength(s2);
 699
 700     count = total = 0;
 701     for (i = 1; i <= shorterSeqLength; i++) // NOTE june29
 702     {
 703         c1 = seqArray[s1][i];
 704         c2 = seqArray[s2][i];
 705         if ((c1 >= 0) && (c1 < userParameters->getMaxAA()))
 706         {
 707             total++;
 708             if (c1 == c2)
 709             {
 710                 count++;
 711             }
 712         }
 713
 714     }
 715
 716     if (total == 0)
 717     {
 718         score = 0;
 719     }
 720     else
 721     {
 722         score = 100.0 *(float)count / (float)total;
 723     }
 724
 725     return (score);
 726 }
 727
 728 void Alignment::debugPrintSequences()
 729 {
 730     cout << std::endl;
 731     for(int i = 0; i < (int)seqArray.size(); i++)
 732     {
 733         for(int j = 0; j < (int)seqArray[i].size(); j++)
 734         {
 735             if(seqArray[i][j] > 9)
 736                 cout << " " << seqArray[i][j];
 737             else
 738                 cout << "  " << seqArray[i][j];
 739         }
 740         cout << std::endl;
 741     }
 742 }
 743
 744 /*
 745  * Note the max_aln_length is now called maxAlignmentLength, and it will be stored
 746  * and calculated in this class. Mostly it is used for allocating arrays. But not always.
 747  */
 748 void Alignment::calculateMaxLengths()
 749 {
 750     maxAlignmentLength = 0;
 751     lengthLongestSequence = 0;
 752     if(seqArray.size() > 0)
 753     {
 754         SeqArray::iterator seqArrayIter = seqArray.begin();
 755
 756         while(seqArrayIter != seqArray.end())
 757         {
 758             // NOTE I needed to change this to >= for a bug I had!!!!!!!
 759             if((int)(*seqArrayIter).size() - 1 >= lengthLongestSequence)
 760             {
 761                 lengthLongestSequence = (*seqArrayIter).size();
 762             }
 763             ++seqArrayIter;
 764         }
 765
 766         if(lengthLongestSequence > 0)
 767         {
 768             maxAlignmentLength = (lengthLongestSequence * 2) - 2;
 769             lengthLongestSequence -= 1; // MADE A CHANGE HERE AS WELL!!
 770         }
 771         else
 772         {
 773             lengthLongestSequence = 0;
 774             maxAlignmentLength = 0;
 775         }
 776
 777     }
 778     else
 779     {
 780         maxAlignmentLength = 0;
 781     }
 782     maxNames = 0;
 783     if(names.size() > 0)
 784     {
 785         vector<string>::iterator nameVecIter = names.begin();
 786
 787         while(nameVecIter != names.end())
 788         {
 789             if((int)(*nameVecIter).size() > maxNames)
 790             {
 791                 maxNames = (*nameVecIter).size();
 792             }
 793             ++nameVecIter;
 794         }
 795         if(maxNames < 10)
 796         {
 797             maxNames = 10;
 798         }
 799     }
 800     else
 801     {
 802         maxNames = 0;
 803     }
 804 }
 805
 806 /*
 807  * This function checks to see if all names are different. It returns true if they
 808  * are all different, and false if there are 2 the same.
 809  */
 810 bool Alignment::checkAllNamesDifferent(string *offendingSeq)
 811 {
 812     bool different = true;
 813     // NOTE I added the + 1 here because, if we had a sequence with a name as blank
 814     // this would be the same as the first one!
 815     vector<string>::iterator namesIter1 = names.begin() + 1;
 816     vector<string>::iterator namesIter2;
 817     int counter1 = 1;
 818     int counter2 = 2;
 819
 820
 821     while(namesIter1 != names.end())
 822     {
 823         namesIter2 = namesIter1 + 1;
 824         while(namesIter2 != names.end())
 825         {
 826             if((*namesIter1).compare(*namesIter2) == 0) // If we have 2 strings the same.
 827             {
 828                 different = false;
 829                 /* 23-03-2007,nige: let someone up the stack deal with this - GUI is too deeply entangled.
 830                  * utilityObject->error("Multiple sequences found with same name '%s' (first %d chars are significant)\n", namesIter1->c_str(), MAXNAMES);
 831                  */
 832
 833                 offendingSeq->assign((*namesIter1));
 834                 clearAlignment();
 835                 return different; // Not all different!
 836             }
 837             namesIter2++;
 838             counter2++;
 839         }
 840         namesIter1++;
 841         counter1++;
 842         counter2 = counter1 + 1;
 843     }
 844     return different;
 845 }
 846
 847
 848
 849
 850
 851
 852 void Alignment::addSecStructMask1(vector<char>* secStructMaskToAdd)
 853 {
 854     secStructMask1 = *secStructMaskToAdd;
 855 }
 856
 857 void Alignment::addSecStructMask2(vector<char>* secStructMaskToAdd)
 858 {
 859     secStructMask2 = *secStructMaskToAdd;
 860 }
 861
 862 void Alignment::addGapPenaltyMask1(vector<char>* gapPenaltyMaskToAdd)
 863 {
 864     gapPenaltyMask1 = *gapPenaltyMaskToAdd;
 865 }
 866
 867 void Alignment::addGapPenaltyMask2(vector<char>* gapPenaltyMaskToAdd)
 868 {
 869     gapPenaltyMask2 = *gapPenaltyMaskToAdd;
 870 }
 871
 872 void Alignment::addSecStructName1(string nameToAdd)
 873 {
 874     secStructName1 = nameToAdd;
 875 }
 876
 877 void Alignment::addSecStructName2(string nameToAdd)
 878 {
 879     secStructName2 = nameToAdd;
 880 }
 881
 882 void Alignment::addSeqWeight(vector<int>* _seqWeight)
 883 {
 884     if(seqWeight.size() == _seqWeight->size())
 885     {
 886         int size = seqWeight.size();
 887
 888         for(int i = 0; i < size; i++)
 889         {
 890             seqWeight[i] = (*_seqWeight)[i];
 891         }
 892     }
 893 }
 894
 895 void Alignment::printSequencesAddedInfo()
 896 {
 897     if(userParameters->getDisplayInfo())
 898     {
 899         int startSeq = userParameters->getProfile2Empty() ? 1:
 900                        profile1NumSeqs + 1;
 901         string dnaFlag = userParameters->getDNAFlag() ? "bp" : "aa";
 902
 903
 904         for(int i = startSeq; i <= numSeqs; i++)
 905         {
 906             cout << "Sequence " << i << ": "
 907                  << std::left << setw(maxNames) << names.at(i)
 908                  << std::right << setw(6) << getSequenceLength(i)
 909                  << " " << dnaFlag << std::endl;
 910         }
 911     }
 912 }
 913
 914 void Alignment::debugPrintOutAlignInfo()
 915 {
 916     for(int i = 1; i <= numSeqs; i++)
 917     {
 918         cout << "seq-no=" << i << ": name="
 919              << std::left << setw(maxNames) << names.at(i)
 920              << " length="
 921              << std::right << setw(6) << getSequenceLength(i)
 922              << std::endl;
 923     }
 924 }
 925
 926 int Alignment::getSequenceLength(int index)
 927 {
 928     return seqArray.at(index).size() - 1;
 929 }
 930
 931 int Alignment::getLengthLongestSequence()
 932 {
 933     int _lengthLongestSequence = 0;
 934
 935     for(int i = 1; i <= numSeqs; i++)
 936     {
 937         if(getSeqLength(i) > _lengthLongestSequence)
 938         {
 939             _lengthLongestSequence = getSeqLength(i);
 940         }
 941     }
 942     return _lengthLongestSequence;
 943 }
 944
 945 /**
 946  * This function is used to find the length of the longest sequence in the range.
 947  */
 948 int Alignment::getLengthLongestSequence(int firstSeq, int lastSeq)
 949 {
 950     int _lengthLongestSequence = 0;
 951
 952     if(firstSeq >= 1 && lastSeq <= numSeqs)
 953     {
 954         for(int i = firstSeq; i <= lastSeq; i++)
 955         {
 956             if(getSeqLength(i) > _lengthLongestSequence)
 957             {
 958                 _lengthLongestSequence = getSeqLength(i);
 959             }
 960         }
 961     }
 962     return _lengthLongestSequence; // Will return 0 if cant check seqs
 963 }
 964
 965 int Alignment::getMaxNames()
 966 {
 967     return maxNames;
 968 }
 969
 970 string Alignment::getSecStructName1()
 971 {
 972     return secStructName1;
 973 }
 974
 975 string Alignment::getSecStructName2()
 976 {
 977     return secStructName2;
 978 }
 979
 980 vector<char>* Alignment::getGapPenaltyMask1()
 981 {
 982     return &gapPenaltyMask1;
 983 }
 984
 985 vector<char>* Alignment::getGapPenaltyMask2()
 986 {
 987     return &gapPenaltyMask2;
 988 }
 989
 990 vector<char>* Alignment::getSecStructMask1()
 991 {
 992     return &secStructMask1;
 993 }
 994
 995 vector<char>* Alignment::getSecStructMask2()
 996 {
 997     return &secStructMask2;
 998 }
 999
1000
1001 int Alignment::getSecStructMask1Element(int index)
1002 {
1003     if(index > 0 && index < (int)secStructMask1.size())
1004     {
1005         return secStructMask1[index];
1006     }
1007     else
1008     {
1009         throw VectorOutOfRange(string("secStructMask1"), index, secStructMask1.size() - 1);
1010     }
1011 }
1012
1013 int Alignment::getSecStructMask2Element(int index)
1014 {
1015     if(index > 0 && index < (int)secStructMask2.size())
1016     {
1017         return secStructMask2[index];
1018     }
1019     else
1020     {
1021         throw VectorOutOfRange(string("secStructMask2"), index, secStructMask2.size() - 1);
1022     }
1023 }
1024
1025 int Alignment::getGapPenaltyMask1Element(int index)
1026 {
1027     if(index > 0 && index < (int)gapPenaltyMask1.size())
1028     {
1029         return gapPenaltyMask1[index];
1030     }
1031     else
1032     {
1033         throw VectorOutOfRange(string("gapPenaltyMask1"), index, gapPenaltyMask1.size() - 1);
1034     }
1035 }
1036
1037 int Alignment::getGapPenaltyMask2Element(int index)
1038 {
1039     if(index > 0 && index < (int)gapPenaltyMask2.size())
1040     {
1041         return gapPenaltyMask2[index];
1042     }
1043     else
1044     {
1045         throw VectorOutOfRange(string("gapPenaltyMask2"), index, gapPenaltyMask2.size() - 1);
1046     }
1047 }
1048
1049 string Alignment::getName(int index)
1050 {
1051     if(index > 0 && index < (int)names.size())
1052     {
1053         return names[index];
1054     }
1055     else
1056     {
1057         throw VectorOutOfRange(string("names"), index, names.size() - 1);
1058     }
1059 }
1060
1061 /**
1062  * Change:
1063  * Mark 25-1-2007: This function was added to allow access to the unique id.
1064  */
1065 unsigned long Alignment::getUniqueId(int seq)
1066 {
1067     if(seq > 0 && seq < (int)sequenceIds.size())
1068     {
1069         return sequenceIds[seq];
1070     }
1071     else
1072     {
1073         throw VectorOutOfRange(string("sequenceIds"), seq, sequenceIds.size() - 1);
1074     }
1075 }
1076
1077 string Alignment::getTitle(int index)
1078 {
1079     if(index > 0 && index < (int)titles.size())
1080     {
1081         return titles[index];
1082     }
1083     else
1084     {
1085         throw VectorOutOfRange(string("titles"), index, titles.size() - 1);
1086     }
1087 }
1088
1089 int Alignment::getOutputIndex(int index)
1090 {
1091     if(index >= 0 && index < (int)outputIndex.size())
1092     {
1093         return outputIndex[index];
1094     }
1095     else
1096     {
1097         throw VectorOutOfRange(string("outputIndex"), index, outputIndex.size() - 1);
1098     }
1099 }
1100
1101 int Alignment::getSeqWeight(int index) const
1102 {
1103     if(index >= 0 && index < (int)seqWeight.size())
1104     {
1105         return seqWeight[index];
1106     }
1107     else
1108     {
1109         throw VectorOutOfRange(string("seqWeight"), index, seqWeight.size() - 1);
1110     }
1111 }
1112
1113
1114 /**
1115  * This function is used by Qt. It is used to calculate the histogram values for the
1116  * widget in clustalQt. The function is in here because it needs access to the SubMatrix
1117  * class.
1118  * @param matNum The number of the matrix to be used in the comparison.
1119  * @return A pointer to a vector containing the histogram values.
1120  */
1121 vector<int>* Alignment::QTcalcHistColumnHeights(int firstSeq, int nSeqs,
1122                                               Array2D<int>* exceptionalRes)
1123 {
1124     int n, i, s, p, r, r1;
1125     int numColumns = getLengthLongestSequence();
1126     int scoreScale = userParameters->getQTScorePlotScale();//5; // From ClustalX.
1127     int scoreCutOff = userParameters->getQTResExceptionCutOff();
1128     bool includeGaps = false;
1129     //short  *mat_xref, *matptr;
1130     float median, mean;
1131     float t, q1, q3, ul;
1132     vector<float> seqdist, sorteddist;
1133     float diff;
1134     vector<int> seqvector;
1135     vector<int> freq;
1136     vector<vector<int> > profile;
1137     int matrix[NUMRES][NUMRES];
1138     int _maxAA = userParameters->getMaxAA();
1139     int _gapPos1 = userParameters->getGapPos1();
1140     histogramColumnHeights.resize(numColumns);
1141     //panel_data data1;
1142     subMatrix->getQTMatrixForHistogram(matrix);
1143
1144     profile.resize(numColumns + 2, vector<int>(_maxAA + 2));
1145     freq.resize(_maxAA + 2);
1146
1147     for(p = 0; p < numColumns; p++)
1148     {
1149         for(r = 0; r < _maxAA; r++)
1150         {
1151             freq[r] = 0;
1152         }
1153         for(s = firstSeq; s < firstSeq + nSeqs; s++)
1154         {
1155             if(p < getSeqLength(s + 1) && seqArray[s + 1][p + 1] >= 0 &&
1156                  seqArray[s + 1][p + 1] < _maxAA)
1157             {
1158                 freq[seqArray[s + 1][p + 1]]++;
1159             }
1160         }
1161         for(r = 0; r < _maxAA; r++)
1162         {
1163             profile[p][r] = 0;
1164             for(r1 = 0; r1 < _maxAA; r1++)
1165             {
1166                 profile[p][r] += freq[r1] * matrix[r1][r];
1167             }
1168             profile[p][r] = static_cast<int>(profile[p][r] /
1169                              static_cast<float>(nSeqs)); // Mark change 17-5-07
1170         }
1171     }
1172
1173     seqvector.resize(_maxAA + 2);
1174     seqdist.resize(nSeqs + 1);
1175     sorteddist.resize(nSeqs + 1);
1176
1177     for(p = 0; p < numColumns; p++)
1178     {
1179         for(s = firstSeq; s < firstSeq + nSeqs; s++)
1180         {
1181             if (p < getSeqLength(s + 1))
1182             {
1183                 for (r = 0; r < _maxAA; r++)
1184                 {
1185                     seqvector[r] = matrix[r][seqArray[s + 1][p + 1]];
1186                 }
1187             }
1188             else
1189             {
1190                 for (r = 0; r < _maxAA; r++)
1191                 {
1192                     seqvector[r] = matrix[r][_gapPos1];
1193                 }
1194             }
1195             seqdist[s - firstSeq] = 0.0;
1196             for(r = 0; r < _maxAA; r++)
1197             {
1198                 diff = profile[p][r] - seqvector[r];
1199                 diff /= 1000.0;
1200                 seqdist[s - firstSeq] += diff * diff;
1201             }
1202             seqdist[s - firstSeq] = sqrt((double)seqdist[s - firstSeq]);
1203         }
1204
1205         // calculate mean,median and rms of seq distances
1206         mean = median = 0.0;
1207         if(includeGaps)
1208         {
1209             for(s = 0; s < nSeqs; s++)
1210             {
1211                 mean += seqdist[s];
1212             }
1213             mean /= nSeqs;
1214             n = nSeqs;
1215             for(s = 0; s < nSeqs; s++)
1216             {
1217                 sorteddist[s] = seqdist[s];
1218             }
1219         }
1220         else
1221         {
1222             n = 0;
1223             for(s = firstSeq; s < firstSeq + nSeqs; s++)
1224             {
1225                 if(p < getSeqLength(s + 1) && seqArray[s + 1][p + 1] >= 0 &&
1226                    seqArray[s + 1][p + 1] < _maxAA)
1227                 {
1228                     mean += seqdist[s - firstSeq];
1229                     n++;
1230                 }
1231             }
1232             if(n > 0)
1233             {
1234                 mean /= n;
1235             }
1236             for(s = firstSeq, i = 0; s < firstSeq + nSeqs; s++)
1237             {
1238                 if(p < getSeqLength(s + 1) && seqArray[s + 1][p + 1] >= 0 &&
1239                    seqArray[s + 1][p + 1] < _maxAA)
1240                 {
1241                     sorteddist[i++] = seqdist[s - firstSeq];
1242                 }
1243             }
1244         }
1245         sortScores(&sorteddist, 0, n - 1);
1246
1247
1248         if(n == 0)
1249         {
1250             median = 0;
1251         }
1252         else if(n % 2 == 0)
1253         {
1254             median = (sorteddist[n / 2 - 1] + sorteddist[n / 2]) / 2.0;
1255         }
1256         else
1257         {
1258             median = sorteddist[n / 2];
1259         }
1260
1261         if(scoreScale <= 5)
1262         {
1263             histogramColumnHeights[p] = static_cast<int>(exp((double)(-mean *
1264                                         (6 - scoreScale) / 4.0)) * 100.0 * n / nSeqs);
1265         }
1266         else
1267         {
1268             histogramColumnHeights[p] = static_cast<int>(exp((double)(-mean /
1269                                         (4.0 * (scoreScale - 4)))) * 100.0 * n / nSeqs);
1270         }
1271
1272         if(n == 0)
1273         {
1274             ul = 0;
1275         }
1276         else
1277         {
1278             t = n/4.0 + 0.5;
1279             if(t - (int)t == 0.5)
1280             {
1281                 q3 = (sorteddist[(int)t] + sorteddist[(int)t + 1]) / 2.0;
1282                 q1 = (sorteddist[n-(int)t] + sorteddist[n - (int)t - 1]) / 2.0;
1283             }
1284             else if(t - (int)t > 0.5)
1285             {
1286                 q3 = sorteddist[(int)t + 1];
1287                 q1 = sorteddist[n - (int)t - 1];
1288             }
1289             else
1290             {
1291                 q3 = sorteddist[(int)t];
1292                 q1 = sorteddist[n - (int)t];
1293             }
1294             if (n < 4)
1295             {
1296                 ul = sorteddist[0];
1297             }
1298             else
1299             {
1300                 ul = q3 + (q3 - q1) * ((float)scoreCutOff / 2.0);
1301             }
1302         }
1303
1304         if((exceptionalRes->getRowSize() >= nSeqs) &&
1305             exceptionalRes->getColSize() >= numColumns)
1306         {
1307             for(int s = firstSeq; s < firstSeq + nSeqs; s++)
1308             {
1309                 if(seqdist[s - firstSeq] > ul && p < getSeqLength(s + 1) &&
1310                    seqArray[s + 1][p + 1] >= 0 &&
1311                    seqArray[s + 1][p + 1] < userParameters->getMaxAA())
1312                 {
1313                     (*exceptionalRes)[s - firstSeq][p] = 1;
1314                 }
1315                 else
1316                 {
1317                     (*exceptionalRes)[s - firstSeq][p] = 0;
1318                 }
1319             }
1320         }
1321     }
1322     return &histogramColumnHeights;
1323 }
1324
1325 void Alignment::sortScores(vector<float>* scores, int f, int l)
1326 {
1327     int i,last;
1328
1329     if(f >= l)
1330         return;
1331
1332     swap(scores, f, (f + l) / 2);
1333     last = f;
1334     for(i = f + 1; i <= l; i++)
1335     {
1336         if((*scores)[i] > (*scores)[f])
1337         {
1338             swap(scores, ++last, i);
1339         }
1340     }
1341     swap(scores, f, last);
1342     sortScores(scores, f, last - 1);
1343     sortScores(scores, last + 1, l);
1344 }
1345
1346 void Alignment::swap(vector<float>* scores, int s1, int s2)
1347 {
1348     float temp;
1349
1350     temp = (*scores)[s1];
1351     (*scores)[s1] = (*scores)[s2];
1352     (*scores)[s2] = temp;
1353 }
1354
1355 int Alignment::searchForString(bool* found, int seq, int beginRes, string search)
1356 {
1357     int lengthSeq = getSeqLength(seq);
1358     if(beginRes > lengthSeq)
1359     {
1360         *found = false;
1361         return beginRes;
1362     }
1363
1364     int res = beginRes;
1365
1366     // First need to convert search into a vector of ints!
1367     vector<int> codedSearch;
1368     int size = search.size();
1369     codedSearch.resize(size);
1370     int code;
1371     for(int i = 0; i < size; i++)
1372     {
1373         code = userParameters->resIndex(userParameters->getAminoAcidCodes(), search[i]);
1374         codedSearch[i] = code;
1375     }
1376
1377     int numSame = 0;
1378     int startPos = -1;
1379     int searchSize = codedSearch.size();
1380     // Now check for the string of ints!!!!
1381     for(; res <= lengthSeq; res++) //- nige: res starts at 1
1382
1383     {
1384         if(seqArray[seq][res] == codedSearch[0])
1385         {
1386             startPos = res; //- nige
1387             for(int i = 0; i < searchSize && (i + res) <= lengthSeq; i++) //- nige: res starts at 1
1388             {
1389                 if(seqArray[seq][res + i] == codedSearch[i])
1390                 {
1391                     numSame++;
1392                 }
1393                 //nige: hack: encoded gap character: see also AlignmentScroll.cpp
1394                 else if (seqArray[seq][res + i] == 31 || seqArray[seq][res + i] == 30)
1395                 {
1396                     res++;
1397                     i--;
1398                 }
1399                 else
1400                 {
1401                     break; // Not the same
1402                 }
1403             }
1404             if(numSame == searchSize)
1405             {
1406                 *found = true;
1407                 return startPos;
1408             }
1409             else
1410             {
1411                 numSame = 0;
1412             }
1413         }
1414     }
1415     *found = false;
1416     return startPos; // Not found!!!!!!!
1417 }
1418
1419 void Alignment::removeGapsFromSelectedSeqs(vector<int>* selected)
1420 {
1421     //getNumSeqs()
1422     int size = selected->size();
1423     int lengthOfSelectedSeq = 0;
1424     int gapPos1 = userParameters->getGapPos1();
1425     int gapPos2 = userParameters->getGapPos2();
1426     int s1;
1427
1428     for(int i = 1; i <= getNumSeqs() && i < size; i++)
1429     {
1430         if((*selected)[i] == 1)
1431         {
1432             // remove gaps from this seq!
1433             lengthOfSelectedSeq = getSeqLength(i);
1434             s1 = 0;
1435             for(int j = 1; j <= lengthOfSelectedSeq; j++)
1436             {
1437                 if(seqArray[i][j] == gapPos1 || seqArray[i][j] == gapPos2)
1438                 {
1439                     continue;
1440                 }
1441                 ++s1;
1442                 seqArray[i][s1] = seqArray[i][j];
1443             }
1444             // Then remove the excess bit at the end of the array
1445             int numExtraElements = lengthOfSelectedSeq - s1;
1446
1447             if((int)seqArray[i].size() > numExtraElements)
1448             {
1449                 for(int k = 0; k < numExtraElements; k++)
1450                 {
1451                     seqArray[i].pop_back();
1452                 }
1453             }
1454         }
1455     }
1456
1457 }
1458
1459 void Alignment::removeGapOnlyColsFromSelectedSeqs(vector<int>* selected)
1460 {
1461     int numGaps = 0;
1462     int NoneSelected = -1;
1463     int numColumns = 0;
1464     int sizeSelected = selected->size();
1465     int firstSeqSelected = NoneSelected;
1466     int gapPos1 = userParameters->getGapPos1();
1467     int gapPos2 = userParameters->getGapPos2();
1468     int k;
1469
1470     for(int i = 1; i < sizeSelected; i++)
1471     {
1472         if((*selected)[i] == 1)
1473         {
1474             numColumns++;
1475             if(firstSeqSelected == -1)
1476             {
1477                 firstSeqSelected = i;
1478             }
1479         }
1480     }
1481
1482     if(firstSeqSelected == NoneSelected)
1483     {
1484         cout << "No Sequences have been selected\n";
1485         return;
1486     }
1487
1488     for(int i = 1; i <= getSeqLength(firstSeqSelected);)
1489     {
1490         numGaps = 0;
1491         for(int j = firstSeqSelected; j < sizeSelected && (*selected)[j] == 1; j++)
1492         {
1493             if(getSeqLength(j) >= i)
1494             {
1495                 if(seqArray[j][i] == gapPos1 || seqArray[j][i] == gapPos2)
1496                 {
1497                     numGaps++;
1498                 }
1499             }
1500         }
1501         if(numGaps == numColumns)
1502         {
1503             //cout << "                removing a gap column\n\n";
1504             for(int j = firstSeqSelected; j < sizeSelected && (*selected)[j] == 1; j++)
1505             {
1506                 for(k = i + 1; k <= getSeqLength(j) + 1 && (int)seqArray[j].size() > k; k++)
1507                 {
1508                     seqArray[j][k - 1] = seqArray[j][k];
1509                 }
1510                 seqArray[j].pop_back(); // Remove the last element!!
1511
1512                 if(getSeqLength(firstSeqSelected) <= 0)
1513                 {
1514                     break;
1515                 }
1516             }
1517         }
1518         else
1519         {
1520             i++;
1521         }
1522     }
1523
1524 }
1525
1526 void Alignment::removeAllGapOnlyColumns(int fSeq, int lSeq, int profileNum)
1527 {
1528     if(fSeq >= lSeq)
1529     {
1530         return;
1531     }
1532     int gapPos1 = userParameters->getGapPos1();
1533     int gapPos2 = userParameters->getGapPos2();
1534
1535     int numGaps = 0;
1536     int numColumns = lSeq - fSeq + 1;
1537     int k;
1538     // We must cheack each column to see if it consists of only '-'
1539     for(int i = 1; i <= getSeqLength(fSeq);)
1540     {
1541         numGaps = 0;
1542         for(int j = fSeq; j <= lSeq; j++)
1543         {
1544             if(getSeqLength(j) >= i)
1545             {
1546                 if(seqArray[j][i] == gapPos1 || seqArray[j][i] == gapPos2)
1547                 {
1548                     numGaps++;
1549                 }
1550             }
1551         }
1552         if(numGaps == numColumns)
1553         {
1554             for(int j = fSeq; j <= lSeq; j++)
1555             {
1556                 for(k = i + 1; k <= getSeqLength(j) + 1; k++)
1557                 {
1558                     seqArray[j][k - 1] = seqArray[j][k];
1559                 }
1560                 seqArray[j].pop_back(); // Remove the last element!!
1561
1562                 if(profileNum == 1)
1563                 {
1564                     int lengthSecStruct = secStructMask1.size();
1565                     int lengthGapMask = gapPenaltyMask1.size();
1566
1567                     for(k = i; k <= getSeqLength(fSeq) && k < lengthSecStruct; k++)
1568                     {
1569                         secStructMask1[k - 1] = secStructMask1[k];
1570                     }
1571                     for(k = i; k <= getSeqLength(fSeq) && k < lengthGapMask; k++)
1572                     {
1573                         gapPenaltyMask1[k - 1] = gapPenaltyMask1[k];
1574                     }
1575                 }
1576
1577                 if(profileNum == 2)
1578                 {
1579                     int lengthSecStruct = secStructMask2.size();
1580                     int lengthGapMask = gapPenaltyMask2.size();
1581
1582                     for(k = i; k <= getSeqLength(fSeq) && k < lengthSecStruct; k++)
1583                     {
1584                         secStructMask2[k - 1] = secStructMask2[k];
1585                     }
1586                     for(k = i; k <= getSeqLength(fSeq) && k < lengthGapMask; k++)
1587                     {
1588                         gapPenaltyMask2[k - 1] = gapPenaltyMask2[k];
1589                     }
1590                 }
1591                 if(getSeqLength(fSeq) <= 0)
1592                 {
1593                     break;
1594                 }
1595             }
1596         }
1597         else
1598         {
1599             i++;
1600         }
1601     }
1602 }
1603
1604 vector<Sequence> Alignment::cutSelectedSequencesFromAlignment(vector<int>* selected)
1605 {
1606     vector<Sequence> cutSequences;
1607     int sizeOfSelected = selected->size();
1608     SeqArray::iterator seqArrayIterator;
1609     vector<string>::iterator namesIterator;
1610     vector<string>::iterator titlesIterator;
1611     vector<int>::iterator seqWeightIterator;
1612     vector<unsigned long>::iterator sequenceIdsIterator;
1613
1614     int newProfile1NumSeqs = profile1NumSeqs;
1615     int profNum = userParameters->getProfileNum();
1616     int prof1NumSeqs = profile1NumSeqs;
1617     int numCutSoFar = 0;
1618     int intialNumSeqs = numSeqs;
1619
1620     for(int i = 1; i < sizeOfSelected && i <= intialNumSeqs; i++)
1621     {
1622         if((*selected)[i] == 1)
1623         {
1624             // Cut the sequence from the alignment!
1625             seqArrayIterator = seqArray.begin() + i - numCutSoFar;
1626             namesIterator = names.begin() + i - numCutSoFar;
1627             titlesIterator = titles.begin() + i - numCutSoFar;
1628             seqWeightIterator = seqWeight.begin() + i - numCutSoFar;
1629             sequenceIdsIterator = sequenceIds.begin() + i - numCutSoFar;
1630             Sequence SeqToCut(&seqArray[i - numCutSoFar], *namesIterator, *titlesIterator,
1631                               *sequenceIdsIterator);
1632
1633             numCutSoFar++;
1634
1635             seqArrayIterator->clear();
1636             seqArray.erase(seqArrayIterator);
1637             names.erase(namesIterator);
1638             titles.erase(titlesIterator);
1639             seqWeight.erase(seqWeightIterator);
1640             sequenceIds.erase(sequenceIdsIterator);
1641
1642             if(numSeqs > 0)
1643             {
1644                 numSeqs--;
1645             }
1646             else
1647             {
1648                 numSeqs = 0;
1649                 break;
1650             }
1651             if(profNum > 0 && i <= prof1NumSeqs)
1652             {
1653                 if(newProfile1NumSeqs > 0)
1654                 {
1655                     newProfile1NumSeqs--;
1656                 }
1657                 else
1658                 {
1659                     newProfile1NumSeqs = 0;
1660                 }
1661             }
1662             cutSequences.push_back(SeqToCut);
1663         }
1664     }
1665     profile1NumSeqs = newProfile1NumSeqs;
1666     setDefaultOutputIndex();
1667     resetAllSeqWeights();
1668     return cutSequences;
1669 }
1670
1671 void Alignment::setDefaultOutputIndex()
1672 {
1673     outputIndex.clear();
1674     outputIndex.resize(numSeqs);
1675     for(int iseq = 1; iseq <= numSeqs; iseq++)
1676     {
1677         outputIndex[iseq - 1] = iseq;
1678     }
1679 }
1680
1681 bool Alignment::removeAllOutsideRange(int beginPos, int endPos)
1682 {
1683     bool ok;
1684     if(beginPos < 0 || endPos > getLengthLongestSequence())
1685     {
1686         return false; // cannot do it!!!!
1687     }
1688
1689     // trim the seqArray
1690     ok = keepPortionOfSeqArray(beginPos, endPos);
1691     if(!ok)
1692     {
1693         cerr << "There was a problem removing a portion of the array\n";
1694         return false;
1695     }
1696
1697     // recalculate the maxLengths
1698     calculateMaxLengths();
1699
1700     // Clear the histogram columns
1701     histogramColumnHeights.clear();
1702
1703     // reset the weights
1704     resetAllSeqWeights();
1705         return true;
1706 }
1707
1708 bool Alignment::keepPortionOfSeqArray(int beginRangeIndex, int endRangeIndex)
1709 {
1710     SeqArray sectionToRealign;
1711     vector<int> emptyVec;
1712     sectionToRealign.push_back(emptyVec); // EMPTY sequence
1713
1714     SeqArray::iterator posToAddTo = sectionToRealign.begin();
1715     // erase from all sequences the range specified here!!!!!
1716     if(beginRangeIndex < 0 || endRangeIndex < 0)
1717     {
1718         return false;
1719     }
1720
1721     SeqArray::iterator mainBeginIt = seqArray.begin() + 1;
1722     SeqArray::iterator mainEndIt = seqArray.end();
1723
1724     vector<int>::iterator begin, end, beginRange, endRange, beginCopyRange, endCopyRange;
1725
1726     for(; mainBeginIt != mainEndIt; mainBeginIt++)
1727     {
1728         vector<int> vecToAdd;
1729         begin = mainBeginIt->begin() + 1;
1730         end = mainBeginIt->end();
1731         beginRange = begin + beginRangeIndex;
1732         endRange = begin + endRangeIndex + 1;
1733         beginCopyRange = beginRange;
1734         endCopyRange = endRange;
1735
1736         // We need to copy all of this into another vector.
1737         if(endCopyRange < end && beginCopyRange < end)
1738         {
1739             vecToAdd.push_back(0);
1740             for(; beginCopyRange != endCopyRange; beginCopyRange++)
1741             {
1742                 vecToAdd.push_back(*beginCopyRange);
1743             }
1744             sectionToRealign.push_back(vecToAdd);
1745         }
1746         else
1747         {
1748             return false;
1749         }
1750
1751         if(endRange < end && beginRange < end)
1752         {
1753             mainBeginIt->erase(beginRange, endRange);
1754         }
1755         else
1756         {
1757             return false;
1758         }
1759     }
1760     clearSeqArray();
1761     seqArray = sectionToRealign;
1762     return true;
1763 }
1764
1765 void Alignment::debugPrintSeqArray(SeqArray* arrayToPrint)
1766 {
1767     // I need to use iterators for everything here.
1768     SeqArray::iterator mainBeginIt = arrayToPrint->begin();
1769     SeqArray::iterator mainEndIt = arrayToPrint->end();
1770     vector<int>::iterator begin, end;
1771     string aaCodes = userParameters->getAminoAcidCodes();
1772
1773     for(; mainBeginIt != mainEndIt; mainBeginIt++)
1774     {
1775         if(mainBeginIt->size() > 0)
1776         {
1777             begin = mainBeginIt->begin() + 1;
1778             end = mainBeginIt->end();
1779             for(; begin != end; begin++)
1780             {
1781                 if(*begin < (int)aaCodes.size())
1782                 {
1783                     cout << aaCodes[*begin];
1784                 }
1785                 else
1786                 {
1787                     cout << "-";
1788                 }
1789             }
1790             cout << "\n";
1791         }
1792     }
1793 }
1794
1795 void Alignment::debugPrintProfile1()
1796 {
1797     cout << "************** PROFILE1 *********************\n";
1798     SeqArray::iterator mainBeginIt = seqArray.begin() + 1;
1799     SeqArray::iterator mainEndIt = mainBeginIt + profile1NumSeqs;
1800     vector<int>::iterator begin, end;
1801     string aaCodes = userParameters->getAminoAcidCodes();
1802
1803     for(; mainBeginIt != mainEndIt; mainBeginIt++)
1804     {
1805         cout << "PROFILE1 SEQ: ";
1806         if(mainBeginIt->size() > 0)
1807         {
1808             begin = mainBeginIt->begin() + 1;
1809             end = mainBeginIt->end();
1810             for(; begin != end; begin++)
1811             {
1812                 if(*begin < (int)aaCodes.size())
1813                 {
1814                     cout << aaCodes[*begin];
1815                 }
1816                 else
1817                 {
1818                     cout << "-";
1819                 }
1820             }
1821             cout << "\n";
1822         }
1823     }
1824 }
1825
1826 void Alignment::debugPrintProfile2()
1827 {
1828     cout << "************** PROFILE2 *********************\n";
1829     SeqArray::iterator mainBeginIt = seqArray.begin() + 1 + profile1NumSeqs;
1830     SeqArray::iterator mainEndIt = seqArray.end();
1831     vector<int>::iterator begin, end;
1832     string aaCodes = userParameters->getAminoAcidCodes();
1833
1834     for(; mainBeginIt != mainEndIt; mainBeginIt++)
1835     {
1836         cout << "PROFILE2 SEQ: ";
1837         if(mainBeginIt->size() > 0)
1838         {
1839             begin = mainBeginIt->begin() + 1;
1840             end = mainBeginIt->end();
1841             for(; begin != end; begin++)
1842             {
1843                 if(*begin < (int)aaCodes.size())
1844                 {
1845                     cout << aaCodes[*begin];
1846                 }
1847                 else
1848                 {
1849                     cout << "-";
1850                 }
1851             }
1852             cout << "\n";
1853         }
1854     }
1855 }
1856
1857 bool Alignment::updateRealignedRange(SeqArray realignedSeqs, int beginPos, int endPos)
1858 {
1859     if(realignedSeqs.size() != seqArray.size())
1860     {
1861         return false;
1862     }
1863     if(beginPos < 0 || endPos < 0)
1864     {
1865         return false;
1866     }
1867
1868     // erase from all sequences the range specified here!!!!!
1869
1870     SeqArray::iterator mainBeginIt = seqArray.begin() + 1;
1871     SeqArray::iterator mainEndIt = seqArray.end();
1872
1873     SeqArray::iterator pasteBeginIt = realignedSeqs.begin() + 1;
1874     SeqArray::iterator pasteEndIt = realignedSeqs.end();
1875
1876     vector<int>::iterator begin, end, beginRange, endRange;
1877
1878     for(; mainBeginIt != mainEndIt && pasteBeginIt != pasteEndIt; mainBeginIt++)
1879     {
1880         vector<int> vecToAdd;
1881         begin = mainBeginIt->begin() + 1;
1882         end = mainBeginIt->end();
1883         beginRange = begin + beginPos;
1884         endRange = begin + endPos + 1;
1885
1886         if(endRange < end && beginRange < end)
1887         {
1888             mainBeginIt->erase(beginRange, endRange);
1889             mainBeginIt->insert(beginRange, pasteBeginIt->begin() + 1,
1890                                 pasteBeginIt->end());
1891         }
1892         else
1893         {
1894             return false;
1895         }
1896         pasteBeginIt++;
1897     }
1898     return true;
1899 }
1900
1901 bool Alignment::reloadAlignment()
1902 {
1903     if(getNumSeqs() <= 0)
1904     {
1905         return false;
1906     }
1907     if(userParameters->getOutputOrder() == INPUT)
1908     {
1909         return true;
1910     }
1911     if((int)outputIndex.size() != getNumSeqs())
1912     {
1913         outputIndex.clear();
1914         return false;
1915     }
1916     vector<int> emptyVec;
1917     string emptyString = "";
1918     SeqArray outputOrderSeqArray;
1919     outputOrderSeqArray.resize(getNumSeqs() + 1);
1920     outputOrderSeqArray[0] = emptyVec;
1921     vector<string> outputOrderNames;
1922     outputOrderNames.resize(getNumSeqs() + 1);
1923     outputOrderNames[0] = emptyString;
1924     vector<string> outputOrderTitles;
1925     outputOrderTitles.resize(getNumSeqs() + 1);
1926     outputOrderTitles[0] = emptyString;
1927     vector<unsigned long> outputOrderSequenceIds;
1928     outputOrderSequenceIds.resize(getNumSeqs() + 1);
1929     outputOrderSequenceIds[0] = 0;
1930
1931     int size = seqArray.size();
1932     if((seqArray.size() != names.size()) || (seqArray.size() != titles.size()) ||
1933         sequenceIds.size() != names.size())
1934     {
1935         return false;
1936     }
1937
1938     int _outIndex;
1939     // Now for each seq,
1940     for(int i = 1; i < size; i++)
1941     {
1942         if(i < (int)outputOrderSeqArray.size() && i - 1 < (int)outputIndex.size() &&
1943            outputIndex[i - 1] < size)
1944         {
1945             _outIndex = outputIndex[i - 1];
1946             outputOrderSeqArray[i] = seqArray[_outIndex];
1947             outputOrderNames[i] = names[_outIndex];
1948             outputOrderTitles[i] = titles[_outIndex];
1949             outputOrderSequenceIds[i] = sequenceIds[_outIndex];
1950         }
1951         else
1952         {
1953             return false;
1954         }
1955     }
1956
1957     // Now we have a copy in the correct order.
1958     // Remove all the elements from the old ones and set them to be these arrays
1959     clearSeqArray();
1960     seqArray = outputOrderSeqArray;
1961     names.clear();
1962     names = outputOrderNames;
1963     titles.clear();
1964     titles = outputOrderTitles;
1965     sequenceIds.clear();
1966     sequenceIds = outputOrderSequenceIds;
1967
1968     return true;
1969 }
1970
1971 const vector<int>* Alignment::getSequenceFromUniqueId(unsigned long id)
1972 {
1973     for(int i = 0; i < (int)sequenceIds.size(); i++)
1974     {
1975         if(sequenceIds[i] == id)
1976         {
1977             return getSequence(i);
1978         }
1979     }
1980
1981     // We have not found it, throw an exception!!!
1982     throw SequenceNotFoundException();
1983 }
1984
1985 /**
1986  * Change:
1987  * Mark 26-1-2007: This function was added to allow access to the unique id.
1988  */
1989 void Alignment::updateSequence(int index, const vector<int>* seq)
1990 {
1991     if(index >= 1 && index < (int)seqArray.size())
1992     {
1993         seqArray[index] = *seq;
1994     }
1995 }
1996
1997 /**
1998  * Change:
1999  * Mark 17-2-2007: This function was added to check if a res is a gap or not.
2000  */
2001 bool Alignment::isGap(int seq, int col) const
2002 {
2003     int res = seqArray[seq][col];
2004     if(res == gapPos1 || res == gapPos2)
2005     {
2006         return true;
2007     }
2008     else
2009     {
2010         return false;
2011     }
2012 }
2013 /**
2014  * This function will be used so that we can create an alignment object from a seqArray.
2015  * This will be used for the tree iteration.
2016  */
2017 void Alignment::addSequences(SeqArray* seqVector)
2018 {
2019     clearAlignment();
2020     numSeqs = seqVector->size() - 1;
2021     vector<int> emptyVec;
2022
2023     seqArray.push_back(emptyVec); // EMPTY sequence
2024     names.push_back(string(""));
2025     titles.push_back(string(""));
2026     sequenceIds.push_back(0);
2027     cout << "\nThere are " << numSeqs << " in the alignment obj\n";
2028     for(int i = 1; i <= numSeqs; i++)
2029     {
2030         ostringstream name;
2031         seqArray.push_back((*seqVector)[i]);
2032         titles.push_back(string(""));
2033         sequenceIds.push_back(utilityObject->getUniqueSequenceIdentifier());
2034         name << "name" << numSeqs;
2035         names.push_back(name.str());
2036     }
2037
2038     calculateMaxLengths();
2039     seqWeight.resize(numSeqs + 1, 100);
2040 }
2041
2042 }
2043