binaries/src/probcons/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
  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
  34   /////////////////////////////////////////////////////////////////
  35   // Sequence::Sequence()
  36   //
  37   // Default constructor.  Does nothing.
  38   /////////////////////////////////////////////////////////////////
  39
  40   Sequence () : isValid (false), header (""), data (NULL), length (0), sequenceLabel (0), inputLabel (0) {}
  41
  42  public:
  43
  44   /////////////////////////////////////////////////////////////////
  45   // Sequence::Sequence()
  46   //
  47   // Constructor.  Reads the sequence from a FileBuffer.
  48   /////////////////////////////////////////////////////////////////
  49
  50   Sequence (FileBuffer &infile, bool stripGaps = false) : isValid (false), header ("~"), data (NULL), length(0), sequenceLabel (0), inputLabel (0) {
  51
  52     // read until the first non-blank line
  53     while (!infile.eof()){
  54       infile.GetLine (header);
  55       if (header.length() != 0) break;
  56     }
  57
  58     // check to make sure that it is a correct header line
  59     if (header[0] == '>'){
  60
  61       // if so, remove the leading ">"
  62       header = header.substr (1);
  63
  64       // remove any leading or trailing white space in the header comment
  65       while (header.length() > 0 && isspace (header[0])) header = header.substr (1);
  66       while (header.length() > 0 && isspace (header[header.length() - 1])) header = header.substr(0, header.length() - 1);
  67
  68       // get ready to read the data[] array; note that data[0] is always '@'
  69       char ch;
  70       data = new SafeVector<char>; assert (data);
  71       data->push_back ('@');
  72
  73       // get a character from the file
  74       while (infile.Get(ch)){
  75
  76         // if we've reached a new comment line, put the character back and stop
  77         if (ch == '>'){ infile.UnGet(); break; }
  78
  79         // skip whitespace
  80         if (isspace (ch)) continue;
  81
  82         // substitute gap character
  83         if (ch == '.') ch = '-';
  84         if (stripGaps && ch == '-') continue;
  85
  86         // check for known characters
  87         if (!((ch >= 'A' && ch <= 'Z') || (ch >= 'a' && ch <= 'z') || ch == '*' || ch == '-')){
  88           cerr << "ERROR: Unknown character encountered: " << ch << endl;
  89           exit (1);
  90         }
  91
  92         // everything's ok so far, so just store this character.
  93         data->push_back(ch);
  94         ++length;
  95       }
  96
  97       // sequence must contain data in order to be valid
  98       isValid = length > 0;
  99       if (!isValid){
 100         delete data;
 101         data = NULL;
 102       }
 103     }
 104   }
 105
 106
 107   /////////////////////////////////////////////////////////////////
 108   // Sequence::Sequence()
 109   //
 110   // Constructor.  Builds a sequence from existing data.  Note
 111   // that the data must use one-based indexing where data[0] should
 112   // be set to '@'.
 113   /////////////////////////////////////////////////////////////////
 114
 115   Sequence (SafeVector<char> *data, string header, int length, int sequenceLabel, int inputLabel) :
 116     isValid (data != NULL), header(header), data(data), length (length), sequenceLabel (sequenceLabel), inputLabel (inputLabel) {
 117       assert (data);
 118       assert ((*data)[0] == '@');
 119   }
 120
 121   /////////////////////////////////////////////////////////////////
 122   // Sequence::Sequence()
 123   //
 124   // Destructor.  Release allocated memory.
 125   /////////////////////////////////////////////////////////////////
 126
 127   ~Sequence (){
 128     if (data){
 129       assert (isValid);
 130       delete data;
 131       data = NULL;
 132       isValid = false;
 133     }
 134   }
 135
 136   /////////////////////////////////////////////////////////////////
 137   // Sequence::GetHeader()
 138   //
 139   // Return the string comment associated with this sequence.
 140   /////////////////////////////////////////////////////////////////
 141
 142   string GetHeader () const {
 143     return header;
 144   }
 145
 146   /////////////////////////////////////////////////////////////////
 147   // Sequence::GetName()
 148   //
 149   // Return the first word of the string comment associated with this sequence.
 150   /////////////////////////////////////////////////////////////////
 151
 152   string GetName () const {
 153     char name[1024];
 154     sscanf (header.c_str(), "%s", name);
 155     return string(name);
 156   }
 157
 158   /////////////////////////////////////////////////////////////////
 159   // Sequence::GetDataPtr()
 160   //
 161   // Return the iterator to data associated with this sequence.
 162   /////////////////////////////////////////////////////////////////
 163
 164   SafeVector<char>::iterator GetDataPtr(){
 165     assert (isValid);
 166     assert (data);
 167     return data->begin();
 168   }
 169
 170   /////////////////////////////////////////////////////////////////
 171   // Sequence::GetPosition()
 172   //
 173   // Return the character at position i.  Recall that the character
 174   // data is stored with one-based indexing.
 175   /////////////////////////////////////////////////////////////////
 176
 177   char GetPosition (int i) const {
 178     assert (isValid);
 179     assert (data);
 180     assert (i >= 1 && i <= length);
 181     return (*data)[i];
 182   }
 183
 184   /////////////////////////////////////////////////////////////////
 185   // Sequence::SetLabel()
 186   //
 187   // Sets the sequence label to i.
 188   /////////////////////////////////////////////////////////////////
 189
 190   void SetLabel (int i){
 191     assert (isValid);
 192     sequenceLabel = i;
 193     inputLabel = i;
 194   }
 195
 196   /////////////////////////////////////////////////////////////////
 197   // Sequence::SetSortLabel()
 198   //
 199   // Sets the sequence sorting label to i.
 200   /////////////////////////////////////////////////////////////////
 201
 202   void SetSortLabel (int i){
 203     assert (isValid);
 204     sequenceLabel = i;
 205   }
 206
 207   /////////////////////////////////////////////////////////////////
 208   // Sequence::GetLabel()
 209   //
 210   // Retrieves the input label.
 211   /////////////////////////////////////////////////////////////////
 212
 213   int GetLabel () const {
 214     assert (isValid);
 215     return inputLabel;
 216   }
 217
 218   /////////////////////////////////////////////////////////////////
 219   // Sequence::GetSortLabel()
 220   //
 221   // Retrieves the sorting label.
 222   /////////////////////////////////////////////////////////////////
 223
 224   int GetSortLabel () const {
 225     assert (isValid);
 226     return sequenceLabel;
 227   }
 228
 229   /////////////////////////////////////////////////////////////////
 230   // Sequence::Fail()
 231   //
 232   // Checks to see if the sequence successfully loaded.
 233   /////////////////////////////////////////////////////////////////
 234
 235   bool Fail () const {
 236     return !isValid;
 237   }
 238
 239   /////////////////////////////////////////////////////////////////
 240   // Sequence::Length()
 241   //
 242   // Returns the length of the sequence.
 243   /////////////////////////////////////////////////////////////////
 244
 245   int GetLength () const {
 246     assert (isValid);
 247     assert (data);
 248     return length;
 249   }
 250
 251   /////////////////////////////////////////////////////////////////
 252   // Sequence::WriteMFA()
 253   //
 254   // Writes the sequence to outfile in MFA format.  Uses numColumns
 255   // columns per line.  If useIndex is set to false, then the
 256   // header is printed as normal, but if useIndex is true, then
 257   // ">S###" is printed where ### represents the sequence label.
 258   /////////////////////////////////////////////////////////////////
 259
 260   void WriteMFA (ostream &outfile, int numColumns, bool useIndex = false) const {
 261     assert (isValid);
 262     assert (data);
 263     assert (!outfile.fail());
 264
 265     // print out heading
 266     if (useIndex)
 267       outfile << ">S" << GetLabel() << endl;
 268     else
 269       outfile << ">" << header << endl;
 270
 271     // print out character data
 272     int ct = 1;
 273     for (; ct <= length; ct++){
 274       outfile << (*data)[ct];
 275       if (ct % numColumns == 0) outfile << endl;
 276     }
 277     if ((ct-1) % numColumns != 0) outfile << endl;
 278   }
 279
 280   /////////////////////////////////////////////////////////////////
 281   // Sequence::Clone()
 282   //
 283   // Returns a new deep copy of the seqeuence.
 284   /////////////////////////////////////////////////////////////////
 285
 286   Sequence *Clone () const {
 287     Sequence *ret = new Sequence();
 288     assert (ret);
 289
 290     ret->isValid = isValid;
 291     ret->header = header;
 292     ret->data = new SafeVector<char>; assert (ret->data);
 293     *(ret->data) = *data;
 294     ret->length = length;
 295     ret->sequenceLabel = sequenceLabel;
 296     ret->inputLabel = inputLabel;
 297
 298     return ret;
 299   }
 300
 301   /////////////////////////////////////////////////////////////////
 302   // Sequence::GetRange()
 303   //
 304   // Returns a new sequence object consisting of a range of
 305   // characters from the current seuquence.
 306   /////////////////////////////////////////////////////////////////
 307
 308   Sequence *GetRange (int start, int end) const {
 309     Sequence *ret = new Sequence();
 310     assert (ret);
 311
 312     assert (start >= 1 && start <= length);
 313     assert (end >= 1 && end <= length);
 314     assert (start <= end);
 315
 316     ret->isValid = isValid;
 317     ret->header = header;
 318     ret->data = new SafeVector<char>; assert (ret->data);
 319     ret->data->push_back ('@');
 320     for (int i = start; i <= end; i++)
 321       ret->data->push_back ((*data)[i]);
 322     ret->length = end - start + 1;
 323     ret->sequenceLabel = sequenceLabel;
 324     ret->inputLabel = inputLabel;
 325
 326     return ret;
 327   }
 328
 329   /////////////////////////////////////////////////////////////////
 330   // Sequence::AddGaps()
 331   //
 332   // Given an SafeVector<char> containing the skeleton for an
 333   // alignment and the identity of the current character, this
 334   // routine will create a new sequence with all necesssary gaps added.
 335   // For instance,
 336   //    alignment = "XXXBBYYYBBYYXX"
 337   //    id = 'X'
 338   // will perform the transformation
 339   //    "ATGCAGTCA" --> "ATGCC---GT--CA"
 340   //                    (XXXBBYYYBBYYXX)
 341   /////////////////////////////////////////////////////////////////
 342
 343   Sequence *AddGaps (SafeVector<char> *alignment, char id){
 344     Sequence *ret = new Sequence();
 345     assert (ret);
 346
 347     ret->isValid = isValid;
 348     ret->header = header;
 349     ret->data = new SafeVector<char>; assert (ret->data);
 350     ret->length = (int) alignment->size();
 351     ret->sequenceLabel = sequenceLabel;
 352     ret->inputLabel = inputLabel;
 353     ret->data->push_back ('@');
 354
 355     SafeVector<char>::iterator dataIter = data->begin() + 1;
 356     for (SafeVector<char>::iterator iter = alignment->begin(); iter != alignment->end(); ++iter){
 357       if (*iter == 'B' || *iter == id){
 358         ret->data->push_back (*dataIter);
 359         ++dataIter;
 360       }
 361       else
 362         ret->data->push_back ('-');
 363     }
 364
 365     return ret;
 366   }
 367
 368   /////////////////////////////////////////////////////////////////
 369   // Sequence::GetString()
 370   //
 371   // Returns the sequence as a string with gaps removed.
 372   /////////////////////////////////////////////////////////////////
 373
 374   string GetString (){
 375     string s = "";
 376     for (int i = 1; i <= length; i++){
 377       if ((*data)[i] != '-') s += (*data)[i];
 378     }
 379     return s;
 380   }
 381
 382
 383   /////////////////////////////////////////////////////////////////
 384   // Sequence::GetMapping()
 385   //
 386   // Returns a SafeVector<int> containing the indices of every
 387   // character in the sequence.  For instance, if the data is
 388   // "ATGCC---GT--CA", the method returns {1,2,3,4,5,9,10,13,14}.
 389   /////////////////////////////////////////////////////////////////
 390
 391   SafeVector<int> *GetMapping () const {
 392     SafeVector<int> *ret = new SafeVector<int>(1, 0);
 393     for (int i = 1; i <= length; i++){
 394       if ((*data)[i] != '-') ret->push_back (i);
 395     }
 396     return ret;
 397   }
 398
 399   /////////////////////////////////////////////////////////////////
 400   // Sequence::Highlight()
 401   //
 402   // Changes all positions with score >= cutoff to upper case and
 403   // all positions with score < cutoff to lower case.
 404   /////////////////////////////////////////////////////////////////
 405
 406   void Highlight (const SafeVector<float> &scores, const float cutoff){
 407     for (int i = 1; i <= length; i++){
 408       if (scores[i-1] >= cutoff)
 409         (*data)[i] = toupper ((*data)[i]);
 410       else
 411         (*data)[i] = tolower ((*data)[i]);
 412     }
 413   }
 414 };
 415
 416 #endif