binaries/src/mafft/extensions/mxscarna_src/probconsRNA/MultiSequence.h

   1 ////////////////////////////////////////////////////////////////
   2 // MultiSequence.h
   3 //
   4 // Utilities for reading/writing multiple sequence data.
   5 /////////////////////////////////////////////////////////////////
   6
   7 #ifndef MULTISEQUENCE_H
   8 #define MULTISEQUENCE_H
   9
  10 #include <cctype>
  11 #include <string>
  12 #include <fstream>
  13 #include <iostream>
  14 #include <sstream>
  15 #include <algorithm>
  16 #include <set>
  17 #include "SafeVector.h"
  18 #include "Sequence.h"
  19 #include "FileBuffer.h"
  20
  21 /////////////////////////////////////////////////////////////////
  22 // MultiSequence
  23 //
  24 // Class for multiple sequence alignment input/output.
  25 /////////////////////////////////////////////////////////////////
  26
  27 namespace MXSCARNA {
  28 class MultiSequence {
  29
  30   SafeVector<Sequence *> *sequences;
  31
  32  public:
  33
  34   /////////////////////////////////////////////////////////////////
  35   // MultiSequence::MultiSequence()
  36   //
  37   // Default constructor.
  38   /////////////////////////////////////////////////////////////////
  39
  40   MultiSequence () : sequences (NULL) {}
  41
  42   /////////////////////////////////////////////////////////////////
  43   // MultiSequence::MultiSequence()
  44   //
  45   // Constructor.  Load MFA from a FileBuffer object.
  46   /////////////////////////////////////////////////////////////////
  47
  48   MultiSequence (FileBuffer &infile) : sequences (NULL) {
  49     LoadMFA (infile);
  50   }
  51
  52   /////////////////////////////////////////////////////////////////
  53   // MultiSequence::MultiSequence()
  54   //
  55   // Constructor.  Load MFA from a filename.
  56   /////////////////////////////////////////////////////////////////
  57
  58   MultiSequence (const string &filename) : sequences (NULL){
  59     LoadMFA (filename);
  60   }
  61
  62   /////////////////////////////////////////////////////////////////
  63   // MultiSequence::~MultiSequence()
  64   //
  65   // Destructor.  Gets rid of sequence objects contained in the
  66   // multiple alignment.
  67   /////////////////////////////////////////////////////////////////
  68
  69   ~MultiSequence(){
  70
  71     // if sequences allocated
  72     if (sequences){
  73
  74       // free all sequences
  75       for (SafeVector<Sequence *>::iterator iter = sequences->begin(); iter != sequences->end(); ++iter){
  76         assert (*iter);
  77         delete *iter;
  78         *iter = NULL;
  79       }
  80
  81       // free sequence vector
  82       delete sequences;
  83       sequences = NULL;
  84     }
  85   }
  86
  87   /////////////////////////////////////////////////////////////////
  88   // MultiSequence::LoadMFA()
  89   //
  90   // Load MFA from a filename.
  91   /////////////////////////////////////////////////////////////////
  92
  93   void LoadMFA (const string &filename, bool stripGaps = false){
  94
  95     // try opening file
  96     FileBuffer infile (filename.c_str());
  97
  98     if (infile.fail()){
  99       cerr << "ERROR: Could not open file '" << filename << "' for reading." << endl;
 100       exit (1);
 101     }
 102
 103     // if successful, then load using other LoadMFA() routine
 104     LoadMFA (infile, stripGaps);
 105
 106     infile.close();
 107   }
 108
 109   /////////////////////////////////////////////////////////////////
 110   // MultiSequence::LoadMFA()
 111   //
 112   // Load MSF from a FileBuffer object.
 113   /////////////////////////////////////////////////////////////////
 114
 115   void ParseMSF (FileBuffer &infile, string header, bool stripGaps = false){
 116
 117     SafeVector<SafeVector<char> *> seqData;
 118     SafeVector<string> seqNames;
 119     SafeVector<int> seqLengths;
 120
 121     istringstream in;
 122     bool valid = true;
 123     bool missingHeader = false;
 124     bool clustalW = false;
 125
 126     // read until data starts
 127     while (!infile.eof() && header.find ("..", 0) == string::npos){
 128       if (header.find ("CLUSTAL", 0) == 0 || header.find ("PROBCONS", 0) == 0){
 129         clustalW = true;
 130         break;
 131       }
 132       infile.GetLine (header);
 133       if (header.find ("//", 0) != string::npos){
 134         missingHeader = true;
 135         break;
 136       }
 137     }
 138
 139     // read until end-of-file
 140     while (valid){
 141       infile.GetLine (header);
 142       if (infile.eof()) break;
 143
 144       string word;
 145       in.clear();
 146       in.str(header);
 147
 148       // check if there's anything on this line
 149       if (in >> word){
 150
 151         // clustalw name parsing
 152         if (clustalW){
 153           if (!isspace(header[0]) && find (seqNames.begin(), seqNames.end(), word) == seqNames.end()){
 154             seqNames.push_back (word);
 155             seqData.push_back (new SafeVector<char>());
 156             seqLengths.push_back (0);
 157             seqData[(int) seqData.size() - 1]->push_back ('@');
 158           }
 159         }
 160
 161         // look for new sequence label
 162         if (word == string ("Name:")){
 163           if (in >> word){
 164             seqNames.push_back (word);
 165             seqData.push_back (new SafeVector<char>());
 166             seqLengths.push_back (0);
 167             seqData[(int) seqData.size() - 1]->push_back ('@');
 168           }
 169           else
 170             valid = false;
 171         }
 172
 173         // check if this is sequence data
 174         else if (find (seqNames.begin(), seqNames.end(), word) != seqNames.end()){
 175           int index = find (seqNames.begin(), seqNames.end(), word) - seqNames.begin();
 176
 177           // read all remaining characters on the line
 178           char ch;
 179           while (in >> ch){
 180             if (isspace (ch)) continue;
 181 //            if (ch >= 'a' && ch <= 'z') ch = ch - 'a' + 'A';
 182             if (ch == '.') ch = '-';
 183             if (stripGaps && ch == '-') continue;
 184 /*
 185             if (!((ch >= 'A' && ch <= 'Z') || ch == '*' || ch == '-')){
 186               cerr << "ERROR: Unknown character encountered: " << ch << endl;
 187               exit (1);
 188             }
 189 */
 190             if (!((ch >= 'A' && ch <= 'Z') || (ch >= 'a' && ch <= 'z') || ch == '*' || ch == '-')){
 191               cerr << "ERROR: Unknown character encountered: " << ch << endl;
 192               exit (1);
 193             }
 194
 195             // everything's ok so far, so just store this character.
 196             seqData[index]->push_back (ch);
 197             seqLengths[index]++;
 198           }
 199         }
 200         else if (missingHeader){
 201           seqNames.push_back (word);
 202           seqData.push_back (new SafeVector<char>());
 203           seqLengths.push_back (0);
 204           seqData[(int) seqData.size() - 1]->push_back ('@');
 205
 206           int index = (int) seqNames.size() - 1;
 207
 208           // read all remaining characters on the line
 209           char ch;
 210           while (in >> ch){
 211             if (isspace (ch)) continue;
 212 //            if (ch >= 'a' && ch <= 'z') ch = ch - 'a' + 'A';
 213             if (ch == '.') ch = '-';
 214             if (stripGaps && ch == '-') continue;
 215
 216             if (!((ch >= 'A' && ch <= 'Z') || ch == '*' || ch == '-')){
 217               cerr << "ERROR: Unknown character encountered: " << ch << endl;
 218               exit (1);
 219             }
 220 /*
 221             if (!((ch >= 'A' && ch <= 'Z') || (ch >= 'a' && ch <= 'z') || ch == '*' || ch == '-')){
 222               cerr << "ERROR: Unknown character encountered: " << ch << endl;
 223               exit (1);
 224             }
 225 */
 226             // everything's ok so far, so just store this character.
 227             seqData[index]->push_back (ch);
 228             seqLengths[index]++;
 229           }
 230         }
 231       }
 232     }
 233
 234     // check for errorsq
 235     if (seqNames.size() == 0){
 236       cerr << "ERROR: No sequences read!" << endl;
 237       exit (1);
 238     }
 239
 240     assert (!sequences);
 241     sequences = new SafeVector<Sequence *>;
 242     for (int i = 0; i < (int) seqNames.size(); i++){
 243       if (seqLengths[i] == 0){
 244         cerr << "ERROR: Sequence of zero length!" << endl;
 245         exit (1);
 246       }
 247       Sequence *seq = new Sequence (seqData[i], seqNames[i], seqLengths[i], i, i);
 248       sequences->push_back (seq);
 249     }
 250   }
 251
 252   /////////////////////////////////////////////////////////////////
 253   // MultiSequence::LoadMFA()
 254   //
 255   // Load MFA from a FileBuffer object.
 256   /////////////////////////////////////////////////////////////////
 257
 258   void LoadMFA (FileBuffer &infile, bool stripGaps = false){
 259
 260     // check to make sure that file reading is ok
 261     if (infile.fail()){
 262       cerr << "ERROR: Error reading file." << endl;
 263       exit (1);
 264     }
 265
 266     // read all sequences
 267     while (true){
 268
 269       // get the sequence label as being the current # of sequences
 270       // NOTE: sequence labels here are zero-based
 271       int index = (!sequences) ? 0 : sequences->size();
 272
 273       // read the sequence
 274       Sequence *seq = new Sequence (infile, stripGaps);
 275       if (seq->Fail()){
 276
 277         // check if alternative file format (i.e. not MFA)
 278         if (index == 0){
 279           string header = seq->GetHeader();
 280           if (header.length() > 0 && header[0] != '>'){
 281             // try MSF format
 282             ParseMSF (infile, header);
 283             break;
 284           }
 285         }
 286
 287         delete seq;
 288         break;
 289       }
 290       seq->SetLabel (index);
 291
 292       // add the sequence to the list of current sequences
 293       if (!sequences) sequences = new SafeVector<Sequence *>;
 294       sequences->push_back (seq);
 295     }
 296
 297     // make sure at least one sequence was read
 298     if (!sequences){
 299       cerr << "ERROR: No sequences read." << endl;
 300       exit (1);
 301     }
 302   }
 303
 304   /////////////////////////////////////////////////////////////////
 305   // MultiSequence::AddSequence()
 306   //
 307   // Add another sequence to an existing sequence list
 308   /////////////////////////////////////////////////////////////////
 309
 310   void AddSequence (Sequence *sequence){
 311     assert (sequence);
 312     assert (!sequence->Fail());
 313
 314     // add sequence
 315     if (!sequences) sequences = new SafeVector<Sequence *>;
 316     sequences->push_back (sequence);
 317   }
 318
 319   /////////////////////////////////////////////////////////////////
 320   // MultiSequence::RemoveSequence()
 321   //
 322   // Remove a sequence from the MultiSequence
 323   /////////////////////////////////////////////////////////////////
 324
 325   void RemoveSequence (int index){
 326     assert (sequences);
 327
 328     assert (index >= 0 && index < (int) sequences->size());
 329     delete (*sequences)[index];
 330
 331     sequences->erase (sequences->begin() + index);
 332   }
 333
 334   /////////////////////////////////////////////////////////////////
 335   // MultiSequence::WriteMFA()
 336   //
 337   // Write MFA to the outfile.  Allows the user to specify the
 338   // number of columns for the output.  Also, useIndices determines
 339   // whether or not the actual sequence comments will be printed
 340   // out or whether the artificially assigned sequence labels will
 341   // be used instead.
 342   /////////////////////////////////////////////////////////////////
 343
 344   void WriteMFA (ostream &outfile, string *ssCons = NULL, int numColumns = 60, bool useIndices = false){
 345     if (!sequences) return;
 346
 347     // loop through all sequences and write them out
 348     for (SafeVector<Sequence *>::iterator iter = sequences->begin(); iter != sequences->end(); ++iter){
 349       (*iter)->WriteMFA (outfile, numColumns, useIndices);
 350     }
 351
 352     int count = 0;
 353     if (ssCons != NULL) {
 354       outfile << ">#=GC SS_cons" << endl;
 355       int length = ssCons->length();
 356       for (int i = 1; i < length; i++ ) {
 357         outfile << ssCons->at(i);
 358         ++count;
 359
 360         if (numColumns <= count) {
 361           outfile << endl;
 362           count = 0;
 363         }
 364
 365       }
 366     }
 367     outfile << endl;
 368   }
 369
 370   void WriteMFAseq (ostream &outfile, int numColumns = 60, bool useIndices = false){
 371     if (!sequences) return;
 372
 373     // loop through all sequences and write them out
 374     for (SafeVector<Sequence *>::iterator iter = sequences->begin(); iter != sequences->end(); ++iter){
 375       (*iter)->WriteMFA (outfile, numColumns, useIndices);
 376     }
 377   }
 378
 379   /////////////////////////////////////////////////////////////////
 380   // MultiSequence::GetAnnotationChar()
 381   //
 382   // Return CLUSTALW annotation for column.
 383   /////////////////////////////////////////////////////////////////
 384
 385   char GetAnnotationChar (SafeVector<char> &column){
 386     SafeVector<int> counts (256, 0);
 387     int allChars = (int) column.size();
 388
 389     for (int i = 0; i < allChars; i++){
 390       counts[(unsigned char) toupper(column[i])]++;
 391     }
 392
 393     allChars -= counts[(unsigned char) '-'];
 394     if (allChars == 1) return ' ';
 395
 396     for (int i = 0; i < 256; i++) if ((char) i != '-' && counts[i] == allChars) return '*';
 397
 398     if (counts[(unsigned char) 'S'] +
 399         counts[(unsigned char) 'T'] +
 400         counts[(unsigned char) 'A'] == allChars)
 401       return ':';
 402
 403     if (counts[(unsigned char) 'N'] +
 404         counts[(unsigned char) 'E'] +
 405         counts[(unsigned char) 'Q'] +
 406         counts[(unsigned char) 'K'] == allChars)
 407       return ':';
 408
 409     if (counts[(unsigned char) 'N'] +
 410         counts[(unsigned char) 'H'] +
 411         counts[(unsigned char) 'Q'] +
 412         counts[(unsigned char) 'K'] == allChars)
 413       return ':';
 414
 415     if (counts[(unsigned char) 'N'] +
 416         counts[(unsigned char) 'D'] +
 417         counts[(unsigned char) 'E'] +
 418         counts[(unsigned char) 'Q'] == allChars)
 419       return ':';
 420
 421     if (counts[(unsigned char) 'Q'] +
 422         counts[(unsigned char) 'H'] +
 423         counts[(unsigned char) 'R'] +
 424         counts[(unsigned char) 'K'] == allChars)
 425       return ':';
 426
 427     if (counts[(unsigned char) 'M'] +
 428         counts[(unsigned char) 'I'] +
 429         counts[(unsigned char) 'L'] +
 430         counts[(unsigned char) 'V'] == allChars)
 431       return ':';
 432
 433     if (counts[(unsigned char) 'M'] +
 434         counts[(unsigned char) 'I'] +
 435         counts[(unsigned char) 'L'] +
 436         counts[(unsigned char) 'F'] == allChars)
 437       return ':';
 438
 439     if (counts[(unsigned char) 'H'] +
 440         counts[(unsigned char) 'Y'] == allChars)
 441       return ':';
 442
 443     if (counts[(unsigned char) 'F'] +
 444         counts[(unsigned char) 'Y'] +
 445         counts[(unsigned char) 'W'] == allChars)
 446       return ':';
 447
 448     if (counts[(unsigned char) 'C'] +
 449         counts[(unsigned char) 'S'] +
 450         counts[(unsigned char) 'A'] == allChars)
 451       return '.';
 452
 453     if (counts[(unsigned char) 'A'] +
 454         counts[(unsigned char) 'T'] +
 455         counts[(unsigned char) 'V'] == allChars)
 456       return '.';
 457
 458     if (counts[(unsigned char) 'S'] +
 459         counts[(unsigned char) 'A'] +
 460         counts[(unsigned char) 'G'] == allChars)
 461       return '.';
 462
 463     if (counts[(unsigned char) 'S'] +
 464         counts[(unsigned char) 'T'] +
 465         counts[(unsigned char) 'N'] +
 466         counts[(unsigned char) 'K'] == allChars)
 467       return '.';
 468
 469     if (counts[(unsigned char) 'S'] +
 470         counts[(unsigned char) 'T'] +
 471         counts[(unsigned char) 'P'] +
 472         counts[(unsigned char) 'A'] == allChars)
 473       return '.';
 474
 475     if (counts[(unsigned char) 'S'] +
 476         counts[(unsigned char) 'G'] +
 477         counts[(unsigned char) 'N'] +
 478         counts[(unsigned char) 'D'] == allChars)
 479       return '.';
 480
 481     if (counts[(unsigned char) 'S'] +
 482         counts[(unsigned char) 'N'] +
 483         counts[(unsigned char) 'D'] +
 484         counts[(unsigned char) 'E'] +
 485         counts[(unsigned char) 'Q'] +
 486         counts[(unsigned char) 'K'] == allChars)
 487       return '.';
 488
 489     if (counts[(unsigned char) 'N'] +
 490         counts[(unsigned char) 'D'] +
 491         counts[(unsigned char) 'E'] +
 492         counts[(unsigned char) 'Q'] +
 493         counts[(unsigned char) 'H'] +
 494         counts[(unsigned char) 'K'] == allChars)
 495       return '.';
 496
 497     if (counts[(unsigned char) 'N'] +
 498         counts[(unsigned char) 'E'] +
 499         counts[(unsigned char) 'H'] +
 500         counts[(unsigned char) 'Q'] +
 501         counts[(unsigned char) 'R'] +
 502         counts[(unsigned char) 'K'] == allChars)
 503       return '.';
 504
 505     if (counts[(unsigned char) 'F'] +
 506         counts[(unsigned char) 'V'] +
 507         counts[(unsigned char) 'L'] +
 508         counts[(unsigned char) 'I'] +
 509         counts[(unsigned char) 'M'] == allChars)
 510       return '.';
 511
 512     if (counts[(unsigned char) 'H'] +
 513         counts[(unsigned char) 'F'] +
 514         counts[(unsigned char) 'Y'] == allChars)
 515       return '.';
 516
 517     return ' ';
 518   }
 519
 520   /////////////////////////////////////////////////////////////////
 521   // MultiSequence::WriteALN()
 522   //
 523   // Write ALN to the outfile.  Allows the user to specify the
 524   // number of columns for the output.
 525   /////////////////////////////////////////////////////////////////
 526
 527   void WriteALN (ostream &outfile, int numColumns = 60){
 528     if (!sequences) return;
 529
 530 //   outfile << "Multplex SCARNA version " << VERSION << " multiple sequence alignment"  << endl;
 531 //   outfile << "PROBCONS version " << VERSION << " multiple sequence alignment" << endl;
 532     outfile << "CLUSTAL W(1.83) multiple sequence alignment" << endl;
 533 //    outfile << "//" << endl;
 534
 535     int longestComment = 0;
 536     SafeVector<SafeVector<char>::iterator> ptrs (GetNumSequences());
 537     SafeVector<int> lengths (GetNumSequences());
 538     for (int i = 0; i < GetNumSequences(); i++){
 539       ptrs[i] = GetSequence (i)->GetDataPtr();
 540       lengths[i] = GetSequence (i)->GetLength();
 541       longestComment = max (longestComment, (int) GetSequence(i)->GetName().length());
 542     }
 543     longestComment += 4;
 544
 545     int writtenChars = 0;
 546     bool allDone = false;
 547
 548     while (!allDone){
 549       outfile << endl;
 550       allDone = true;
 551
 552       // loop through all sequences and write them out
 553       for (int i = 0; i < GetNumSequences(); i++){
 554
 555         if (writtenChars < lengths[i]){
 556           outfile << GetSequence(i)->GetName();
 557           for (int j = 0; j < longestComment - (int) GetSequence(i)->GetName().length(); j++)
 558             outfile << ' ';
 559
 560           for (int j = 0; j < numColumns; j++){
 561             if (writtenChars + j < lengths[i])
 562               outfile << ptrs[i][writtenChars + j + 1];
 563             else
 564               break;
 565           }
 566
 567           outfile << endl;
 568
 569           if (writtenChars + numColumns < lengths[i]) allDone = false;
 570         }
 571       }
 572
 573       // write annotation line
 574       for (int j = 0; j < longestComment; j++)
 575         outfile << ' ';
 576
 577       for (int j = 0; j < numColumns; j++){
 578         SafeVector<char> column;
 579
 580         for (int i = 0; i < GetNumSequences(); i++)
 581           if (writtenChars + j < lengths[i])
 582             column.push_back (ptrs[i][writtenChars + j + 1]);
 583
 584         if (column.size() > 0)
 585           outfile << GetAnnotationChar (column);
 586       }
 587
 588       outfile << endl;
 589       writtenChars += numColumns;
 590     }
 591     outfile << endl;
 592   }
 593
 594   ////////////////////////////////////////////////////////////////
 595   // MultiSequence::WriteWEB();
 596   //
 597   // Write ALN to the outfile.  Allows the user to specify the
 598   // number of columns for the output.
 599   ///////////////////////////////////////////////////////////////
 600    void WriteWEB (ostream &outfile, string *ssCons = NULL, int numColumns = 60, bool useIndices = false){
 601     if (!sequences) return;
 602
 603     // loop through all sequences and write them out
 604     for (SafeVector<Sequence *>::iterator iter = sequences->begin(); iter != sequences->end(); ++iter){
 605       (*iter)->WriteWEB (outfile, numColumns, useIndices);
 606     }
 607
 608     // write conservation
 609     outfile << "<conservation>" << endl;
 610         int longestComment = 0;
 611     SafeVector<SafeVector<char>::iterator> ptrs (GetNumSequences());
 612     SafeVector<int> lengths (GetNumSequences());
 613     for (int i = 0; i < GetNumSequences(); i++){
 614       ptrs[i] = GetSequence (i)->GetDataPtr();
 615       lengths[i] = GetSequence (i)->GetLength();
 616       longestComment = max (longestComment, (int) GetSequence(i)->GetName().length());
 617     }
 618     longestComment += 4;
 619
 620     int writtenChars = 0;
 621     bool allDone = false;
 622
 623     while (!allDone){
 624 //      outfile << endl;
 625       allDone = true;
 626
 627       // loop through all sequences and write them out
 628       for (int i = 0; i < GetNumSequences(); i++){
 629
 630         if (writtenChars < lengths[i]){
 631 //        outfile << GetSequence(i)->GetName();
 632           for (int j = 0; j < longestComment - (int) GetSequence(i)->GetName().length(); j++)
 633 //          outfile << ' ';
 634
 635           for (int j = 0; j < numColumns; j++){
 636               if (writtenChars + j < lengths[i]);
 637 //            outfile << ptrs[i][writtenChars + j + 1];
 638             else
 639               break;
 640           }
 641
 642 //        outfile << endl;
 643
 644           if (writtenChars + numColumns < lengths[i]) allDone = false;
 645         }
 646       }
 647
 648       // write annotation line
 649 //      for (int j = 0; j < longestComment; j++)
 650 //      outfile << ' ';
 651
 652       for (int j = 0; j < numColumns; j++){
 653         SafeVector<char> column;
 654
 655         for (int i = 0; i < GetNumSequences(); i++)
 656           if (writtenChars + j < lengths[i])
 657             column.push_back (ptrs[i][writtenChars + j + 1]);
 658
 659         if (column.size() > 0)
 660           outfile << GetAnnotationChar (column);
 661       }
 662
 663 //      outfile << endl;
 664       writtenChars += numColumns;
 665     }
 666     outfile << endl;
 667     outfile << "</conservation>" << endl;
 668
 669     // write structure information
 670     if (ssCons != NULL) {
 671       outfile << "<structure>" << endl;
 672       int length = ssCons->length();
 673       for (int i = 1; i < length; i++ ) {
 674         outfile << ssCons->at(i);
 675       }
 676       outfile << endl;
 677       outfile << "</structure>" << endl;
 678
 679       // add coordinate information 06/09/14
 680       outfile << "<coordinate>" << endl;
 681
 682       int segmentPos = 1;
 683       for (int i = 1; i < length; i++) {
 684         int count = 0;
 685
 686         if ( ssCons->at(i) == '(' ) {
 687           ++count;
 688
 689           int j = i;
 690           while (count != 0) {
 691             char ch = ssCons->at(++j);
 692             if      (ch == '(')
 693               ++count;
 694             else if (ch == ')')
 695               --count;
 696           }
 697
 698           outfile << "<segment position=\"" << segmentPos++ << "\" starts=\""
 699                   << i << "\"" << " ends=\"" << j << "\"/>" << endl;
 700
 701         }
 702       }
 703     }
 704     outfile << "</coordinate>" << endl;
 705
 706     outfile << "<mxscarna>" << endl;
 707     WriteMXSCARNA (outfile, ssCons);
 708     outfile << "</mxscarna>" << endl;
 709
 710     outfile << "<aln>" << endl;
 711     WriteALN (outfile);
 712     outfile << "</aln>" << endl;
 713
 714     outfile << "<mfa>" << endl;
 715     WriteMFA (outfile, ssCons);
 716     outfile << "</mfa>" << endl;
 717
 718     outfile << "<stockholm>" << endl;
 719     WriteWebSTOCKHOLM (outfile, ssCons);
 720     outfile << "</stockholm>" << endl;
 721   }
 722
 723   ////////////////////////////////////////////////////////////////
 724   // MultiSequence::WriteSTOCKHOLM();
 725   //
 726   // Write STOCKHOLM to the outfile.  Allows the user to specify the
 727   // number of columns for the output.
 728   ///////////////////////////////////////////////////////////////
 729   void WriteSTOCKHOLM (ostream &outfile, string *ssCons = NULL, int numColumns = 60) {
 730     if (!sequences) return;
 731
 732         outfile << "# STOCKHOLM 1.0" << endl;
 733
 734     int longestComment = 0;
 735     SafeVector<SafeVector<char>::iterator> ptrs (GetNumSequences());
 736     SafeVector<int> lengths (GetNumSequences());
 737     for (int i = 0; i < GetNumSequences(); i++){
 738       ptrs[i] = GetSequence (i)->GetDataPtr();
 739       lengths[i] = GetSequence (i)->GetLength();
 740       longestComment = max (longestComment, (int) GetSequence(i)->GetName().length());
 741     }
 742     longestComment += 4;
 743
 744     int writtenChars = 0;
 745     bool allDone = false;
 746
 747     while (!allDone){
 748       outfile << endl;
 749       allDone = true;
 750
 751       // loop through all sequences and write them out
 752       for (int i = 0; i < GetNumSequences(); i++){
 753
 754         if (writtenChars < lengths[i]){
 755           outfile << GetSequence(i)->GetName();
 756           for (int j = 0; j < longestComment - (int) GetSequence(i)->GetName().length(); j++)
 757             outfile << ' ';
 758
 759           for (int j = 0; j < numColumns; j++){
 760             if (writtenChars + j < lengths[i])
 761                 if (ptrs[i][writtenChars + j + 1] != '-')
 762                   outfile << ptrs[i][writtenChars + j + 1];
 763                 else
 764                   outfile << ".";
 765             else
 766               break;
 767           }
 768
 769           outfile << endl;
 770
 771           if (writtenChars + numColumns < lengths[i]) allDone = false;
 772         }
 773       }
 774
 775       // write ssCons
 776
 777       if (ssCons != NULL) {
 778           outfile << "#=GC SS_cons";
 779           int lengthSScons = 12;
 780           for (int j = 0; j < longestComment - lengthSScons; j++)
 781               outfile << ' ';
 782
 783           for (int j = 0; j < numColumns; j++) {
 784               if (ssCons->at(writtenChars + j + 1) == '(')
 785                 outfile << "<";
 786               else if (ssCons->at(writtenChars + j + 1) == ')')
 787                 outfile << ">";
 788               else
 789                 outfile << ".";
 790               if ((unsigned int)writtenChars + j + 1 >= ssCons->length() - 1)
 791                   break;
 792           }
 793           outfile << endl;
 794       }
 795
 796       writtenChars += numColumns;
 797     }
 798     outfile << "//";
 799     outfile << endl;
 800   }
 801
 802     ////////////////////////////////////////////////////////////////
 803   // MultiSequence::WriteSTOCKHOLM();
 804   //
 805   // Write STOCKHOLM to the outfile.  Allows the user to specify the
 806   // number of columns for the output.
 807   ///////////////////////////////////////////////////////////////
 808   void WriteWebSTOCKHOLM (ostream &outfile, string *ssCons = NULL, int numColumns = 60) {
 809     if (!sequences) return;
 810
 811         outfile << "# STOCKHOLM 1.0" << endl;
 812
 813     int longestComment = 0;
 814     SafeVector<SafeVector<char>::iterator> ptrs (GetNumSequences());
 815     SafeVector<int> lengths (GetNumSequences());
 816     for (int i = 0; i < GetNumSequences(); i++){
 817       ptrs[i] = GetSequence (i)->GetDataPtr();
 818       lengths[i] = GetSequence (i)->GetLength();
 819       longestComment = max (longestComment, (int) GetSequence(i)->GetName().length());
 820     }
 821     longestComment += 4;
 822
 823     int writtenChars = 0;
 824     bool allDone = false;
 825
 826     while (!allDone){
 827       outfile << endl;
 828       allDone = true;
 829
 830       // loop through all sequences and write them out
 831       for (int i = 0; i < GetNumSequences(); i++){
 832
 833         if (writtenChars < lengths[i]){
 834           outfile << GetSequence(i)->GetName();
 835           for (int j = 0; j < longestComment - (int) GetSequence(i)->GetName().length(); j++)
 836             outfile << ' ';
 837
 838           for (int j = 0; j < numColumns; j++){
 839             if (writtenChars + j < lengths[i])
 840                 if (ptrs[i][writtenChars + j + 1] != '-')
 841                   outfile << ptrs[i][writtenChars + j + 1];
 842                 else
 843                   outfile << ".";
 844             else
 845               break;
 846           }
 847
 848           outfile << endl;
 849
 850           if (writtenChars + numColumns < lengths[i]) allDone = false;
 851         }
 852       }
 853
 854       // write ssCons
 855
 856       if (ssCons != NULL) {
 857           outfile << "#=GC SS_cons";
 858           int lengthSScons = 12;
 859           for (int j = 0; j < longestComment - lengthSScons; j++)
 860               outfile << ' ';
 861
 862           for (int j = 0; j < numColumns; j++) {
 863               outfile << ssCons->at(writtenChars + j + 1);
 864
 865               if ((unsigned int)writtenChars + j + 1 >= ssCons->length() - 1)
 866                   break;
 867           }
 868           outfile << endl;
 869       }
 870
 871       writtenChars += numColumns;
 872     }
 873     outfile << "//";
 874     outfile << endl;
 875   }
 876
 877   ////////////////////////////////////////////////////////////////
 878   // MultiSequence::WriteMXSCARNA();
 879   //
 880   // Write MXSCARNA to the outfile.  Allows the user to specify the
 881   // number of columns for the output.
 882   ///////////////////////////////////////////////////////////////
 883   void WriteMXSCARNA (ostream &outfile, string *ssCons = NULL, int numColumns = 60){
 884     if (!sequences) return;
 885
 886    outfile << "Multplex SCARNA version " << VERSION << " multiple sequence alignment"  << endl;
 887
 888     int longestComment = 0;
 889     SafeVector<SafeVector<char>::iterator> ptrs (GetNumSequences());
 890     SafeVector<int> lengths (GetNumSequences());
 891     for (int i = 0; i < GetNumSequences(); i++){
 892       ptrs[i] = GetSequence (i)->GetDataPtr();
 893       lengths[i] = GetSequence (i)->GetLength();
 894       longestComment = max (longestComment, (int) GetSequence(i)->GetName().length());
 895     }
 896     longestComment += 4;
 897
 898     int writtenChars = 0;
 899     bool allDone = false;
 900
 901     while (!allDone){
 902       outfile << endl;
 903       allDone = true;
 904
 905       // loop through all sequences and write them out
 906       for (int i = 0; i < GetNumSequences(); i++){
 907
 908         if (writtenChars < lengths[i]){
 909           outfile << GetSequence(i)->GetName();
 910           for (int j = 0; j < longestComment - (int) GetSequence(i)->GetName().length(); j++)
 911             outfile << ' ';
 912
 913           for (int j = 0; j < numColumns; j++){
 914             if (writtenChars + j < lengths[i])
 915               outfile << ptrs[i][writtenChars + j + 1];
 916             else
 917               break;
 918           }
 919
 920           outfile << endl;
 921
 922           if (writtenChars + numColumns < lengths[i]) allDone = false;
 923         }
 924       }
 925
 926       // write ssCons
 927       if (ssCons != NULL) {
 928           outfile << "ss_cons";
 929           int lengthSScons = 7;
 930           for (int j = 0; j < longestComment - lengthSScons; j++)
 931               outfile << ' ';
 932
 933           for (int j = 0; j < numColumns; j++) {
 934               outfile << ssCons->at(writtenChars + j + 1);
 935               if ((unsigned int)writtenChars + j + 1 >= ssCons->length() - 1)
 936                   break;
 937           }
 938           outfile << endl;
 939       }
 940
 941       // write annotation line
 942       for (int j = 0; j < longestComment; j++)
 943         outfile << ' ';
 944
 945       for (int j = 0; j < numColumns; j++){
 946         SafeVector<char> column;
 947
 948         for (int i = 0; i < GetNumSequences(); i++)
 949           if (writtenChars + j < lengths[i])
 950             column.push_back (ptrs[i][writtenChars + j + 1]);
 951
 952         if (column.size() > 0)
 953           outfile << GetAnnotationChar (column);
 954       }
 955
 956       outfile << endl;
 957       writtenChars += numColumns;
 958     }
 959     outfile << endl;
 960   }
 961
 962   /////////////////////////////////////////////////////////////////
 963   // MultiSequence::GetSequence()
 964   //
 965   // Retrieve a sequence from the MultiSequence object.
 966   /////////////////////////////////////////////////////////////////
 967
 968   Sequence* GetSequence (int i){
 969     assert (sequences);
 970     assert (0 <= i && i < (int) sequences->size());
 971
 972     return (*sequences)[i];
 973   }
 974
 975   /////////////////////////////////////////////////////////////////
 976   // MultiSequence::GetSequence()
 977   //
 978   // Retrieve a sequence from the MultiSequence object
 979   // (const version).
 980   /////////////////////////////////////////////////////////////////
 981
 982   const Sequence* GetSequence (int i) const {
 983     assert (sequences);
 984     assert (0 <= i && i < (int) sequences->size());
 985
 986     return (*sequences)[i];
 987   }
 988
 989   /////////////////////////////////////////////////////////////////
 990   // MultiSequence::GetNumSequences()
 991   //
 992   // Returns the number of sequences in the MultiSequence.
 993   /////////////////////////////////////////////////////////////////
 994
 995   int GetNumSequences () const {
 996     if (!sequences) return 0;
 997     return (int) sequences->size();
 998   }
 999
1000   /////////////////////////////////////////////////////////////////
1001   // MultiSequence::SortByHeader()
1002   //
1003   // Organizes the sequences according to their sequence headers
1004   // in ascending order.
1005   /////////////////////////////////////////////////////////////////
1006
1007   void SortByHeader () {
1008     assert (sequences);
1009
1010     // a quick and easy O(n^2) sort
1011     for (int i = 0; i < (int) sequences->size()-1; i++){
1012       for (int j = i+1; j < (int) sequences->size(); j++){
1013         if ((*sequences)[i]->GetHeader() > (*sequences)[j]->GetHeader())
1014           swap ((*sequences)[i], (*sequences)[j]);
1015       }
1016     }
1017   }
1018
1019   /////////////////////////////////////////////////////////////////
1020   // MultiSequence::SortByLabel()
1021   //
1022   // Organizes the sequences according to their sequence labels
1023   // in ascending order.
1024   /////////////////////////////////////////////////////////////////
1025
1026   void SortByLabel () {
1027     assert (sequences);
1028
1029     // a quick and easy O(n^2) sort
1030     for (int i = 0; i < (int) sequences->size()-1; i++){
1031       for (int j = i+1; j < (int) sequences->size(); j++){
1032         if ((*sequences)[i]->GetSortLabel() > (*sequences)[j]->GetSortLabel())
1033           swap ((*sequences)[i], (*sequences)[j]);
1034       }
1035     }
1036   }
1037
1038   /////////////////////////////////////////////////////////////////
1039   // MultiSequence::SaveOrdering()
1040   //
1041   // Relabels sequences so as to preserve the current ordering.
1042   /////////////////////////////////////////////////////////////////
1043
1044   void SaveOrdering () {
1045     assert (sequences);
1046
1047     for (int i = 0; i < (int) sequences->size(); i++)
1048       (*sequences)[i]->SetSortLabel (i);
1049   }
1050
1051   /////////////////////////////////////////////////////////////////
1052   // MultiSequence::Project()
1053   //
1054   // Given a set of indices, extract all sequences from the current
1055   // MultiSequence object whose index is included in the set.
1056   // Then, project the multiple alignments down to the desired
1057   // subset, and return the projection as a new MultiSequence
1058   // object.
1059   /////////////////////////////////////////////////////////////////
1060
1061   MultiSequence *Project (const set<int> &indices){
1062     SafeVector<SafeVector<char>::iterator> oldPtrs (indices.size());
1063     SafeVector<SafeVector<char> *> newPtrs (indices.size());
1064
1065     assert (indices.size() != 0);
1066
1067     // grab old data
1068     int i = 0;
1069     for (set<int>::const_iterator iter = indices.begin(); iter != indices.end(); ++iter){
1070       oldPtrs[i++] = GetSequence (*iter)->GetDataPtr();
1071     }
1072
1073     // compute new length
1074     int oldLength = GetSequence (*indices.begin())->GetLength();
1075     int newLength = 0;
1076     for (i = 1; i <= oldLength; i++){
1077
1078       // check to see if there is a gap in every sequence of the set
1079       bool found = false;
1080       for (int j = 0; !found && j < (int) indices.size(); j++)
1081         found = (oldPtrs[j][i] != '-');
1082
1083       // if not, then this column counts towards the sequence length
1084       if (found) newLength++;
1085     }
1086
1087     // build new alignments
1088     for (i = 0; i < (int) indices.size(); i++){
1089       newPtrs[i] = new SafeVector<char>(); assert (newPtrs[i]);
1090       newPtrs[i]->push_back ('@');
1091     }
1092
1093     // add all needed columns
1094     for (i = 1; i <= oldLength; i++){
1095
1096       // make sure column is not gapped in all sequences in the set
1097       bool found = false;
1098       for (int j = 0; !found && j < (int) indices.size(); j++)
1099         found = (oldPtrs[j][i] != '-');
1100
1101       // if not, then add it
1102       if (found){
1103         for (int j = 0; j < (int) indices.size(); j++)
1104           newPtrs[j]->push_back (oldPtrs[j][i]);
1105       }
1106     }
1107
1108     // wrap sequences in MultiSequence object
1109     MultiSequence *ret = new MultiSequence();
1110     i = 0;
1111     for (set<int>::const_iterator iter = indices.begin(); iter != indices.end(); ++iter){
1112       ret->AddSequence (new Sequence (newPtrs[i++], GetSequence (*iter)->GetHeader(), newLength,
1113                                       GetSequence (*iter)->GetSortLabel(), GetSequence (*iter)->GetLabel()));
1114     }
1115
1116     return ret;
1117   }
1118 };
1119 }
1120 #endif