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

   1 /////////////////////////////////////////////////////////////////
   2 // Sequence.h
   3 //
   4 // Class for reading/manipulating single sequence character data.
   5 /////////////////////////////////////////////////////////////////
   6
   7 #ifndef __SEQUENCE_H__
   8 #define __SEQUENCE_H__
   9
  10 #include <string>
  11 #include <fstream>
  12 #include <iostream>
  13 #include <cctype>
  14 #include <cstdlib>
  15 #include "SafeVector.h"
  16 #include "FileBuffer.h"
  17
  18 /////////////////////////////////////////////////////////////////
  19 // Sequence
  20 //
  21 // Class for storing sequence information.
  22 /////////////////////////////////////////////////////////////////
  23 namespace MXSCARNA {
  24 class Sequence {
  25
  26   bool isValid;                // a boolean indicating whether the sequence data is valid or not
  27   string header;               // string containing the comment line of the FASTA file
  28   SafeVector<char> *data;      // pointer to character data
  29   int length;                  // length of the sequence
  30   int sequenceLabel;           // integer sequence label, typically to indicate the ordering of sequences
  31                                //   in a Multi-FASTA file
  32   int inputLabel;              // position of sequence in original input
  33   float weight;
  34
  35   /////////////////////////////////////////////////////////////////
  36   // Sequence::Sequence()
  37   //
  38   // Default constructor.  Does nothing.
  39   /////////////////////////////////////////////////////////////////
  40
  41   Sequence () : isValid (false), header (""), data (NULL), length (0), sequenceLabel (0), inputLabel (0) {}
  42
  43  public:
  44
  45   /////////////////////////////////////////////////////////////////
  46   // Sequence::Sequence()
  47   //
  48   // Constructor.  Reads the sequence from a FileBuffer.
  49   /////////////////////////////////////////////////////////////////
  50
  51   Sequence (FileBuffer &infile, bool stripGaps = false) : isValid (false), header ("~"), data (NULL), length(0), sequenceLabel (0), inputLabel (0) {
  52
  53     // read until the first non-blank line
  54     while (!infile.eof()){
  55       infile.GetLine (header);
  56       if (header.length() != 0) break;
  57     }
  58
  59     // check to make sure that it is a correct header line
  60     if (header[0] == '>'){
  61
  62       // if so, remove the leading ">"
  63       header = header.substr (1);
  64
  65       // remove any leading or trailing white space in the header comment
  66       while (header.length() > 0 && isspace (header[0])) header = header.substr (1);
  67       while (header.length() > 0 && isspace (header[header.length() - 1])) header = header.substr(0, header.length() - 1);
  68
  69       // get ready to read the data[] array; note that data[0] is always '@'
  70       char ch;
  71       data = new SafeVector<char>; assert (data);
  72       data->push_back ('@');
  73
  74       // get a character from the file
  75       while (infile.Get(ch)){
  76
  77         // if we've reached a new comment line, put the character back and stop
  78         if (ch == '>'){ infile.UnGet(); break; }
  79
  80         // skip whitespace
  81         if (isspace (ch)) continue;
  82
  83         // substitute gap character
  84         if (ch == '.') ch = '-';
  85         if (stripGaps && ch == '-') continue;
  86
  87         // check for known characters
  88         if (!((ch >= 'A' && ch <= 'Z') || (ch >= 'a' && ch <= 'z') || ch == '*' || ch == '-')){
  89           cerr << "ERROR: Unknown character encountered: " << ch << endl;
  90           exit (1);
  91         }
  92
  93         // everything's ok so far, so just store this character.
  94         data->push_back(ch);
  95         ++length;
  96       }
  97
  98       // sequence must contain data in order to be valid
  99       isValid = length > 0;
 100       if (!isValid){
 101         delete data;
 102         data = NULL;
 103       }
 104     }
 105   }
 106
 107
 108   /////////////////////////////////////////////////////////////////
 109   // Sequence::Sequence()
 110   //
 111   // Constructor.  Builds a sequence from existing data.  Note
 112   // that the data must use one-based indexing where data[0] should
 113   // be set to '@'.
 114   /////////////////////////////////////////////////////////////////
 115
 116   Sequence (SafeVector<char> *data, string header, int length, int sequenceLabel, int inputLabel) :
 117     isValid (data != NULL), header(header), data(data), length (length), sequenceLabel (sequenceLabel), inputLabel (inputLabel) {
 118       assert (data);
 119       assert ((*data)[0] == '@');
 120   }
 121
 122   /////////////////////////////////////////////////////////////////
 123   // Sequence::Sequence()
 124   //
 125   // Destructor.  Release allocated memory.
 126   /////////////////////////////////////////////////////////////////
 127
 128   ~Sequence (){
 129     if (data){
 130       assert (isValid);
 131       delete data;
 132       data = NULL;
 133       isValid = false;
 134     }
 135   }
 136
 137     void SetWeight(float myWeight) {
 138       weight = myWeight;
 139   }
 140   float GetWeight() const {
 141       return weight;
 142   }
 143
 144   /////////////////////////////////////////////////////////////////
 145   // Sequence::GetHeader()
 146   //
 147   // Return the string comment associated with this sequence.
 148   /////////////////////////////////////////////////////////////////
 149
 150   string GetHeader () const {
 151     return header;
 152   }
 153
 154   /////////////////////////////////////////////////////////////////
 155   // Sequence::GetName()
 156   //
 157   // Return the first word of the string comment associated with this sequence.
 158   /////////////////////////////////////////////////////////////////
 159
 160   string GetName () const {
 161     char name[1024];
 162     sscanf (header.c_str(), "%s", name);
 163     return string(name);
 164   }
 165
 166   /////////////////////////////////////////////////////////////////
 167   // Sequence::GetDataPtr()
 168   //
 169   // Return the iterator to data associated with this sequence.
 170   /////////////////////////////////////////////////////////////////
 171
 172   SafeVector<char>::iterator GetDataPtr(){
 173     assert (isValid);
 174     assert (data);
 175     return data->begin();
 176   }
 177
 178   /////////////////////////////////////////////////////////////////
 179   // Sequence::GetPosition()
 180   //
 181   // Return the character at position i.  Recall that the character
 182   // data is stored with one-based indexing.
 183   /////////////////////////////////////////////////////////////////
 184
 185   char GetPosition (int i) const {
 186     assert (isValid);
 187     assert (data);
 188     assert (i >= 0 && i <= length);
 189     return (*data)[i];
 190   }
 191
 192   /////////////////////////////////////////////////////////////////
 193   // Sequence::SetLabel()
 194   //
 195   // Sets the sequence label to i.
 196   /////////////////////////////////////////////////////////////////
 197
 198   void SetLabel (int i){
 199     assert (isValid);
 200     sequenceLabel = i;
 201     inputLabel = i;
 202   }
 203
 204   /////////////////////////////////////////////////////////////////
 205   // Sequence::SetSortLabel()
 206   //
 207   // Sets the sequence sorting label to i.
 208   /////////////////////////////////////////////////////////////////
 209
 210   void SetSortLabel (int i){
 211     assert (isValid);
 212     sequenceLabel = i;
 213   }
 214
 215   /////////////////////////////////////////////////////////////////
 216   // Sequence::GetLabel()
 217   //
 218   // Retrieves the input label.
 219   /////////////////////////////////////////////////////////////////
 220
 221   int GetLabel () const {
 222     assert (isValid);
 223     return inputLabel;
 224   }
 225
 226   /////////////////////////////////////////////////////////////////
 227   // Sequence::GetSortLabel()
 228   //
 229   // Retrieves the sorting label.
 230   /////////////////////////////////////////////////////////////////
 231
 232   int GetSortLabel () const {
 233     assert (isValid);
 234     return sequenceLabel;
 235   }
 236
 237   /////////////////////////////////////////////////////////////////
 238   // Sequence::Fail()
 239   //
 240   // Checks to see if the sequence successfully loaded.
 241   /////////////////////////////////////////////////////////////////
 242
 243   bool Fail () const {
 244     return !isValid;
 245   }
 246
 247   /////////////////////////////////////////////////////////////////
 248   // Sequence::Length()
 249   //
 250   // Returns the length of the sequence.
 251   /////////////////////////////////////////////////////////////////
 252
 253   int GetLength () const {
 254     assert (isValid);
 255     assert (data);
 256     return length;
 257   }
 258
 259   /////////////////////////////////////////////////////////////////
 260   // Sequence::WriteMFA()
 261   //
 262   // Writes the sequence to outfile in MFA format.  Uses numColumns
 263   // columns per line.  If useIndex is set to false, then the
 264   // header is printed as normal, but if useIndex is true, then
 265   // ">S###" is printed where ### represents the sequence label.
 266   /////////////////////////////////////////////////////////////////
 267
 268   void WriteMFA (ostream &outfile, int numColumns, bool useIndex = false) const {
 269     assert (isValid);
 270     assert (data);
 271     assert (!outfile.fail());
 272
 273     // print out heading
 274     if (useIndex)
 275       outfile << ">S" << GetLabel() << endl;
 276     else
 277       outfile << ">" << header << endl;
 278
 279     // print out character data
 280     int ct = 1;
 281     for (; ct <= length; ct++){
 282       outfile << (*data)[ct];
 283       if (ct % numColumns == 0) outfile << endl;
 284     }
 285     if ((ct-1) % numColumns != 0) outfile << endl;
 286   }
 287
 288   /////////////////////////////////////////////////////////////////
 289   // Sequence::WriteWEB()
 290   //
 291   // output for web interfase based on Sequence::WriteMFA()
 292   /////////////////////////////////////////////////////////////////
 293
 294   void WriteWEB (ostream &outfile, int numColumns, bool useIndex = false) const {
 295     assert (isValid);
 296     assert (data);
 297     assert (!outfile.fail());
 298
 299     outfile << "<php ref=\"" << GetLabel() << "\">" << endl;
 300     outfile << "<name>" << endl;
 301     // print out heading
 302     if (useIndex)
 303       outfile << "S" << GetLabel() << endl;
 304     else
 305       outfile << "" << header << endl;
 306
 307     outfile << "</name>" << endl;
 308
 309     // print out character data
 310     outfile << "<sequence>" << endl;
 311     int ct = 1;
 312     for (; ct <= length; ct++){
 313       outfile << (*data)[ct];
 314       if (ct % numColumns == 0) outfile << endl;
 315     }
 316     if ((ct-1) % numColumns != 0) outfile << endl;
 317
 318     outfile << "</sequence>" << endl;
 319     outfile << "</php>" << endl;
 320   }
 321
 322   /////////////////////////////////////////////////////////////////
 323   // Sequence::Clone()
 324   //
 325   // Returns a new deep copy of the seqeuence.
 326   /////////////////////////////////////////////////////////////////
 327
 328   Sequence *Clone () const {
 329     Sequence *ret = new Sequence();
 330     assert (ret);
 331
 332     ret->isValid = isValid;
 333     ret->header = header;
 334     ret->data = new SafeVector<char>; assert (ret->data);
 335     *(ret->data) = *data;
 336     ret->length = length;
 337     ret->sequenceLabel = sequenceLabel;
 338     ret->inputLabel = inputLabel;
 339     ret->weight = weight;
 340
 341     return ret;
 342   }
 343
 344   /////////////////////////////////////////////////////////////////
 345   // Sequence::GetRange()
 346   //
 347   // Returns a new sequence object consisting of a range of
 348   // characters from the current seuquence.
 349   /////////////////////////////////////////////////////////////////
 350
 351   Sequence *GetRange (int start, int end) const {
 352     Sequence *ret = new Sequence();
 353     assert (ret);
 354
 355     assert (start >= 1 && start <= length);
 356     assert (end >= 1 && end <= length);
 357     assert (start <= end);
 358
 359     ret->isValid = isValid;
 360     ret->header = header;
 361     ret->data = new SafeVector<char>; assert (ret->data);
 362     ret->data->push_back ('@');
 363     for (int i = start; i <= end; i++)
 364       ret->data->push_back ((*data)[i]);
 365     ret->length = end - start + 1;
 366     ret->sequenceLabel = sequenceLabel;
 367     ret->inputLabel = inputLabel;
 368
 369     return ret;
 370   }
 371
 372   /////////////////////////////////////////////////////////////////
 373   // Sequence::AddGaps()
 374   //
 375   // Given an SafeVector<char> containing the skeleton for an
 376   // alignment and the identity of the current character, this
 377   // routine will create a new sequence with all necesssary gaps added.
 378   // For instance,
 379   //    alignment = "XXXBBYYYBBYYXX"
 380   //    id = 'X'
 381   // will perform the transformation
 382   //    "ATGCAGTCA" --> "ATGCC---GT--CA"
 383   //                    (XXXBBYYYBBYYXX)
 384   /////////////////////////////////////////////////////////////////
 385
 386   Sequence *AddGaps (SafeVector<char> *alignment, char id){
 387     Sequence *ret = new Sequence();
 388     assert (ret);
 389
 390     ret->isValid = isValid;
 391     ret->header = header;
 392     ret->data = new SafeVector<char>; assert (ret->data);
 393     ret->length = (int) alignment->size();
 394     ret->sequenceLabel = sequenceLabel;
 395     ret->inputLabel = inputLabel;
 396     ret->data->push_back ('@');
 397
 398     SafeVector<char>::iterator dataIter = data->begin() + 1;
 399     for (SafeVector<char>::iterator iter = alignment->begin(); iter != alignment->end(); ++iter){
 400       if (*iter == 'B' || *iter == id){
 401         ret->data->push_back (*dataIter);
 402         ++dataIter;
 403       }
 404       else
 405         ret->data->push_back ('-');
 406     }
 407
 408     return ret;
 409   }
 410
 411   /////////////////////////////////////////////////////////////////
 412   // Sequence::AddGaps()
 413   //
 414   // Given an SafeVector<char> containing the skeleton for an
 415   // alignment and the identity of the current character, this
 416   // routine will create a new sequence with all necesssary gaps added.
 417   // For instance,
 418   //    alignment = "XXXBBYYYBBYYXX"
 419   //    id = 'X'
 420   // will perform the transformation
 421   //    "ATGCAGTCA" --> "ATGCC---GT--CA"
 422   //                    (XXXBBYYYBBYYXX)
 423   /////////////////////////////////////////////////////////////////
 424   Sequence *AddGapsReverse (SafeVector<char> *alignment, char id){
 425     Sequence *ret = new Sequence();
 426     assert (ret);
 427
 428     ret->isValid = isValid;
 429     ret->header = header;
 430     ret->data = new SafeVector<char>; assert (ret->data);
 431     ret->length = (int) alignment->size();
 432     ret->sequenceLabel = sequenceLabel;
 433     ret->inputLabel = inputLabel;
 434     ret->data->push_back ('@');
 435
 436     SafeVector<char>::iterator dataIter = data->begin() + 1;
 437     for (SafeVector<char>::reverse_iterator iter = alignment->rbegin(); iter != alignment->rend(); ++iter){
 438       if (*iter == 'B' || *iter == id){
 439         ret->data->push_back (*dataIter);
 440         ++dataIter;
 441       }
 442       else
 443         ret->data->push_back ('-');
 444     }
 445
 446     return ret;
 447   }
 448
 449
 450   /////////////////////////////////////////////////////////////////
 451   // Sequence::GetString()
 452   //
 453   // Returns the sequence as a string with gaps removed.
 454   /////////////////////////////////////////////////////////////////
 455
 456   string GetString (){
 457     string s = " ";
 458     for (int i = 1; i <= length; i++){
 459       if ((*data)[i] != '-') s += (*data)[i];
 460     }
 461     return s;
 462   }
 463
 464
 465   /////////////////////////////////////////////////////////////////
 466   // Sequence::GetMapping()
 467   //
 468   // Returns a SafeVector<int> containing the indices of every
 469   // character in the sequence.  For instance, if the data is
 470   // "ATGCC---GT--CA", the method returns {1,2,3,4,5,9,10,13,14}.
 471   /////////////////////////////////////////////////////////////////
 472
 473   SafeVector<int> *GetMapping () const {
 474     SafeVector<int> *ret = new SafeVector<int>(1, 0);
 475     for (int i = 1; i <= length; i++){
 476       if ((*data)[i] != '-') ret->push_back (i);
 477     }
 478     return ret;
 479   }
 480
 481   /////////////////////////////////////////////////////////////////
 482   // Sequence::GetMappingNumber()
 483   //
 484   // Returns a SafeVector<int> containing the indices of every
 485   // character in the sequence.  For instance, if the data is
 486   // "ATGCC---GT--CA", the method returns {1,2,3,4,5,0,0,0,6,7,0,0,8,9}.
 487   /////////////////////////////////////////////////////////////////
 488   SafeVector<int> *GetMappingNumber () const {
 489       SafeVector<int> *ret = new SafeVector<int>(1, 0);
 490       int count = 0;
 491       for(int i = 1; i <= length; i++) {
 492           if((*data)[i] != '-') ret->push_back(++count);
 493           else                  ret->push_back(0);
 494       }
 495       return ret;
 496   }
 497
 498   /////////////////////////////////////////////////////////////////
 499   // Sequence::Highlight()
 500   //
 501   // Changes all positions with score >= cutoff to upper case and
 502   // all positions with score < cutoff to lower case.
 503   /////////////////////////////////////////////////////////////////
 504
 505   void Highlight (const SafeVector<float> &scores, const float cutoff){
 506     for (int i = 1; i <= length; i++){
 507       if (scores[i-1] >= cutoff)
 508         (*data)[i] = toupper ((*data)[i]);
 509       else
 510         (*data)[i] = tolower ((*data)[i]);
 511     }
 512   }
 513 };
 514 }
 515 #endif // __SQUENCE_HPP__