src/utils/bamtools_fasta.cpp

   1 // ***************************************************************************
   2 // bamtools_fasta.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 13 July 2010
   7 // ---------------------------------------------------------------------------
   8 // Provides FASTA reading/indexing functionality.
   9 // ***************************************************************************
  10
  11 #include <utils/bamtools_fasta.h>
  12 using namespace BamTools;
  13
  14 #include <cstdio>
  15 #include <cstdlib>
  16 #include <cstring>
  17 #include <fstream>
  18 #include <iostream>
  19 #include <sstream>
  20 #include <vector>
  21 using namespace std;
  22
  23 #ifdef _MSC_VER
  24         #define ftello _ftelli64
  25         #define fseeko _fseeki64
  26 #endif
  27
  28 struct Fasta::FastaPrivate {
  29
  30     struct FastaIndexData {
  31         string  Name;
  32         int32_t Length;
  33         int64_t Offset;
  34         int32_t LineLength;
  35         int32_t ByteLength; // LineLength + newline character(s) - varies on OS where file was generated
  36     };
  37
  38     // data members
  39     FILE* Stream;
  40     bool IsOpen;
  41
  42     FILE* IndexStream;
  43     bool HasIndex;
  44     bool IsIndexOpen;
  45
  46     vector<FastaIndexData> Index;
  47
  48     // ctor
  49     FastaPrivate(void);
  50     ~FastaPrivate(void);
  51
  52     // 'public' API methods
  53     bool Close(void);
  54     bool CreateIndex(const string& indexFilename);
  55     bool GetBase(const int& refId, const int& position, char& base);
  56     bool GetSequence(const int& refId, const int& start, const int& stop, string& sequence);
  57     bool Open(const string& filename, const string& indexFilename);
  58
  59     // internal methods
  60     private:
  61         void Chomp(char* sequence);
  62         bool GetNameFromHeader(const string& header, string& name);
  63         bool GetNextHeader(string& header);
  64         bool GetNextSequence(string& sequence);
  65         bool LoadIndexData(void);
  66         bool Rewind(void);
  67         bool WriteIndexData(void);
  68 };
  69
  70 Fasta::FastaPrivate::FastaPrivate(void)
  71     : IsOpen(false)
  72     , HasIndex(false)
  73     , IsIndexOpen(false)
  74 { }
  75
  76 Fasta::FastaPrivate::~FastaPrivate(void) {
  77     Close();
  78 }
  79
  80 // remove any trailing newlines
  81 void Fasta::FastaPrivate::Chomp(char* sequence) {
  82
  83     static const int CHAR_LF = 10;
  84     static const int CHAR_CR = 13;
  85
  86     size_t seqLength = strlen(sequence);
  87     if ( seqLength == 0 ) return;
  88     --seqLength; // ignore null terminator
  89
  90     while ( sequence[seqLength] == CHAR_LF ||
  91             sequence[seqLength] == CHAR_CR
  92           )
  93     {
  94         sequence[seqLength] = 0;
  95         --seqLength;
  96         if (seqLength < 0)
  97             break;
  98     }
  99 }
 100
 101 bool Fasta::FastaPrivate::Close(void) {
 102
 103     // close fasta file
 104     if ( IsOpen ) {
 105         fclose(Stream);
 106         IsOpen = false;
 107     }
 108
 109     // close index file
 110     if ( HasIndex && IsIndexOpen ) {
 111         fclose(IndexStream);
 112         HasIndex = false;
 113         IsIndexOpen = false;
 114     }
 115
 116     // return success
 117     return true;
 118 }
 119
 120 bool Fasta::FastaPrivate::CreateIndex(const string& indexFilename) {
 121
 122     // check that file is open
 123     if ( !IsOpen ) {
 124         cerr << "FASTA error : cannot create index, FASTA file not open" << endl;
 125         return false;
 126     }
 127
 128     // rewind FASTA file
 129     if ( !Rewind() ) {
 130         cerr << "FASTA error : could not rewind FASTA file" << endl;
 131         return false;
 132     }
 133
 134     // clear out prior index data
 135     Index.clear();
 136
 137     // -------------------------------------------
 138     // calculate lineLength & byteLength
 139
 140     int lineLength = 0;
 141     int byteLength = 0;
 142
 143     // skip over header
 144     char buffer[1024];
 145     if ( fgets(buffer, 1024, Stream) == 0 ) {
 146         cerr << "FASTA error : could not read from file" << endl;
 147         return false;
 148     }
 149     if ( feof(Stream) ) return false;
 150     if ( buffer[0] != '>' ) {
 151         cerr << "FASTA error : expected header ('>'), instead : " << buffer[0] << endl;
 152         return false;
 153     }
 154
 155     // read in first line of sequence
 156     char c = fgetc(Stream);
 157     while ( (c >= 0) && (c != '\n') ) {
 158         ++byteLength;
 159         if (isgraph(c)) ++lineLength;
 160         c = fgetc(Stream);
 161     }
 162     ++byteLength; // store newline
 163
 164     // rewind FASTA file
 165     if ( !Rewind() ) {
 166         cerr << "FASTA error : could not rewind FASTA file" << endl;
 167         return false;
 168     }
 169
 170     // iterate through fasta entries
 171     int currentId   = 0;
 172     string header   = "";
 173     string sequence = "";
 174     while ( GetNextHeader(header) ) {
 175
 176         // ---------------------------
 177         // build index entry data
 178         FastaIndexData data;
 179
 180         // store file offset of beginning of DNA sequence (after header)
 181         data.Offset = ftello(Stream);
 182
 183         // parse header, store sequence name in data.Name
 184         if ( !GetNameFromHeader(header, data.Name) ) {
 185             cerr << "FASTA error : could not parse read name from FASTA header" << endl;
 186             return false;
 187         }
 188
 189         // retrieve FASTA sequence
 190         if ( !GetNextSequence(sequence) ) {
 191             cerr << "FASTA error : could not read in next sequence from FASTA file" << endl;
 192             return false;
 193         }
 194
 195         // store sequence length & line/byte lengths
 196         data.Length = sequence.length();
 197         data.LineLength = lineLength;
 198         data.ByteLength = byteLength;
 199
 200         // store index entry
 201         Index.push_back(data);
 202
 203         // update ref Id
 204         ++currentId;
 205     }
 206
 207     // open index file
 208     if ( !indexFilename.empty() ) {
 209         IndexStream = fopen(indexFilename.c_str(), "wb");
 210         if ( !IndexStream ) {
 211             cerr << "FASTA error : Could not open " << indexFilename << " for writing." << endl;
 212             return false;
 213         }
 214         IsIndexOpen = true;
 215     }
 216
 217     // write index data
 218     if ( !WriteIndexData() ) return false;
 219     HasIndex = true;
 220
 221     // close index file
 222     fclose(IndexStream);
 223     IsIndexOpen = false;
 224
 225     // return succes status
 226     return true;
 227 }
 228
 229 bool Fasta::FastaPrivate::GetBase(const int& refId, const int& position, char& base) {
 230
 231     // make sure FASTA file is open
 232     if ( !IsOpen ) {
 233         cerr << "FASTA error : file not open for reading" << endl;
 234         return false;
 235     }
 236
 237     // use index if available
 238     if ( HasIndex && !Index.empty() ) {
 239
 240         // validate reference id
 241         if ( (refId < 0) || (refId >= (int)Index.size()) ) {
 242             cerr << "FASTA error: invalid refId specified: " << refId << endl;
 243             return false;
 244         }
 245
 246         // retrieve reference index data
 247         const FastaIndexData& referenceData = Index.at(refId);
 248
 249         // validate position
 250         if ( (position < 0) || (position > referenceData.Length) ) {
 251             cerr << "FASTA error: invalid position specified: " << position << endl;
 252             return false;
 253         }
 254
 255         // seek to beginning of sequence data
 256         if ( fseeko(Stream, referenceData.Offset, SEEK_SET) != 0 ) {
 257             cerr << "FASTA error : could not sek in file" << endl;
 258             return false;
 259         }
 260
 261         // retrieve sequence
 262         string sequence = "";
 263         if ( !GetNextSequence(sequence) ) {
 264             cerr << "FASTA error : could not retrieve base from FASTA file" << endl;
 265             return false;
 266         }
 267
 268         // set base & return success
 269         base = sequence.at(position);
 270         return true;
 271     }
 272
 273     // else plow through sequentially
 274     else {
 275
 276         // rewind FASTA file
 277         if ( !Rewind() ) {
 278             cerr << "FASTA error : could not rewind FASTA file" << endl;
 279             return false;
 280         }
 281
 282         // iterate through fasta entries
 283         int currentId = 0;
 284         string header = "";
 285         string sequence = "";
 286
 287         // get first entry
 288         GetNextHeader(header);
 289         GetNextSequence(sequence);
 290
 291         while ( currentId != refId ) {
 292             GetNextHeader(header);
 293             GetNextSequence(sequence);
 294             ++currentId;
 295         }
 296
 297         // get desired base from sequence
 298         // TODO: error reporting on invalid position
 299         if ( currentId == refId && (sequence.length() >= (size_t)position) ) {
 300             base = sequence.at(position);
 301             return true;
 302         }
 303
 304         // could not get sequence
 305         return false;
 306     }
 307
 308     // return success
 309     return true;
 310 }
 311
 312 bool Fasta::FastaPrivate::GetNameFromHeader(const string& header, string& name) {
 313
 314     // get rid of the leading greater than sign
 315     string s = header.substr(1);
 316
 317     // extract the first non-whitespace segment
 318     char* pName = (char*)s.data();
 319     unsigned int nameLen = (unsigned int)s.size();
 320
 321     unsigned int start = 0;
 322     while ( (pName[start] == 32) || (pName[start] == 9) || (pName[start] == 10) || (pName[start] == 13) ) {
 323         start++;
 324         if ( start == nameLen )
 325             break;
 326     }
 327
 328     unsigned int stop  = start;
 329     if ( stop < nameLen ) {
 330         while( (pName[stop] != 32) && (pName[stop] != 9) && (pName[stop] != 10) && (pName[stop] != 13) ) {
 331             stop++;
 332             if ( stop == nameLen )
 333                 break;
 334         }
 335     }
 336
 337     if ( start == stop ) {
 338         cerr << "FASTA error : could not parse read name from FASTA header" << endl;
 339         return false;
 340     }
 341
 342     name = s.substr(start, stop - start).c_str();
 343     return true;
 344 }
 345
 346 bool Fasta::FastaPrivate::GetNextHeader(string& header) {
 347
 348     // validate input stream
 349     if ( !IsOpen || feof(Stream) )
 350         return false;
 351
 352     // read in header line
 353     char buffer[1024];
 354     if ( fgets(buffer, 1024, Stream) == 0 ) {
 355         cerr << "FASTA error : could not read from file" << endl;
 356         return false;
 357     }
 358
 359     // make sure it's a FASTA header
 360     if ( buffer[0] != '>' ) {
 361         cerr << "FASTA error : expected header ('>'), instead : " << buffer[0] << endl;
 362         return false;
 363     }
 364
 365     // import buffer contents to header string
 366     stringstream headerBuffer("");
 367     headerBuffer << buffer;
 368     header = headerBuffer.str();
 369
 370     // return success
 371     return true;
 372 }
 373
 374 bool Fasta::FastaPrivate::GetNextSequence(string& sequence) {
 375
 376     // validate input stream
 377     if ( !IsOpen || feof(Stream) )
 378         return false;
 379
 380     // read in sequence
 381     char buffer[1024];
 382     ostringstream seqBuffer("");
 383     while(true) {
 384
 385         char ch = fgetc(Stream);
 386         ungetc(ch, Stream);
 387         if( (ch == '>') || feof(Stream) )
 388               break;
 389
 390         if ( fgets(buffer, 1024, Stream) == 0 ) {
 391             cerr << "FASTA error : could not read from file" << endl;
 392             return false;
 393         }
 394
 395         Chomp(buffer);
 396         seqBuffer << buffer;
 397     }
 398
 399     // import buffer contents to sequence string
 400     sequence = seqBuffer.str();
 401
 402     // return success
 403     return true;
 404 }
 405
 406 bool Fasta::FastaPrivate::GetSequence(const int& refId, const int& start, const int& stop, string& sequence) {
 407
 408     // make sure FASTA file is open
 409     if ( !IsOpen ) {
 410         cerr << "FASTA error : file not open for reading" << endl;
 411         return false;
 412     }
 413
 414     // use index if available
 415     if ( HasIndex && !Index.empty() ) {
 416
 417         // validate reference id
 418         if ( (refId < 0) || (refId >= (int)Index.size()) ) {
 419             cerr << "FASTA error: invalid refId specified: " << refId << endl;
 420             return false;
 421         }
 422
 423         // retrieve reference index data
 424         const FastaIndexData& referenceData = Index.at(refId);
 425
 426         // validate stop position
 427         if ( (start < 0) || (start > stop) || (stop > referenceData.Length) ) {
 428             cerr << "FASTA error: invalid start/stop positions specified: " << start << ", " << stop << endl;
 429             return false;
 430         }
 431
 432         // seek to beginning of sequence data
 433         if ( fseeko(Stream, referenceData.Offset, SEEK_SET) != 0 ) {
 434             cerr << "FASTA error : could not sek in file" << endl;
 435             return false;
 436         }
 437
 438         // retrieve full sequence
 439         string fullSequence = "";
 440         if ( !GetNextSequence(fullSequence) ) {
 441             cerr << "FASTA error : could not retrieve sequence from FASTA file" << endl;
 442             return false;
 443         }
 444
 445         // set sub-sequence & return success
 446         const int seqLength = (stop - start) + 1;
 447         sequence = fullSequence.substr(start, seqLength);
 448         return true;
 449     }
 450
 451     // else plow through sequentially
 452     else {
 453
 454         // rewind FASTA file
 455         if ( !Rewind() ) {
 456             cerr << "FASTA error : could not rewind FASTA file" << endl;
 457             return false;
 458         }
 459
 460         // iterate through fasta entries
 461         int currentId = 0;
 462         string header = "";
 463         string fullSequence = "";
 464
 465         // get first entry
 466         GetNextHeader(header);
 467         GetNextSequence(fullSequence);
 468
 469         while ( currentId != refId ) {
 470             GetNextHeader(header);
 471             GetNextSequence(fullSequence);
 472             ++currentId;
 473         }
 474
 475         // get desired substring from sequence
 476         // TODO: error reporting on invalid start/stop positions
 477         if ( currentId == refId && (fullSequence.length() >= (size_t)stop) ) {
 478             const int seqLength = (stop - start) + 1;
 479             sequence = fullSequence.substr(start, seqLength);
 480             return true;
 481         }
 482
 483         // could not get sequence
 484         return false;
 485     }
 486
 487     // return success
 488     return true;
 489 }
 490
 491 bool Fasta::FastaPrivate::LoadIndexData(void) {
 492
 493     // skip if no index file available
 494     if ( !IsIndexOpen ) return false;
 495
 496     // clear any prior index data
 497     Index.clear();
 498
 499     char buffer[1024];
 500     stringstream indexBuffer;
 501     while ( true ) {
 502
 503         char c = fgetc(IndexStream);
 504         if ( (c == '\n') || feof(IndexStream) ) break;
 505         ungetc(c, IndexStream);
 506
 507         // clear index buffer
 508         indexBuffer.str("");
 509
 510         // read line from index file
 511         if ( fgets(buffer, 1024, IndexStream) == 0 ) {
 512             cerr << "FASTA LoadIndexData() error : could not read from index file" << endl;
 513             HasIndex = false;
 514             return false;
 515         }
 516
 517         // store line in indexBuffer
 518         indexBuffer << buffer;
 519
 520         // retrieve fasta index data from line
 521         FastaIndexData data;
 522         indexBuffer >> data.Name;
 523         indexBuffer >> data.Length;
 524         indexBuffer >> data.Offset;
 525         indexBuffer >> data.LineLength;
 526         indexBuffer >> data.ByteLength;
 527
 528         // store index entry
 529         Index.push_back(data);
 530     }
 531
 532     return true;
 533 }
 534
 535 bool Fasta::FastaPrivate::Open(const string& filename, const string& indexFilename) {
 536
 537     bool success = true;
 538
 539     // open FASTA filename
 540     Stream = fopen(filename.c_str(), "rb");
 541     if ( !Stream ) {
 542         cerr << "FASTA error: Could not open " << filename << " for reading" << endl;
 543         return false;
 544     }
 545     IsOpen = true;
 546     success &= IsOpen;
 547
 548     // open index file if it exists
 549     if ( !indexFilename.empty() ) {
 550         IndexStream = fopen(indexFilename.c_str(), "rb");
 551         if ( !IndexStream ) {
 552             cerr << "FASTA error : Could not open " << indexFilename << " for reading." << endl;
 553             return false;
 554         }
 555         IsIndexOpen = true;
 556         success &= IsIndexOpen;
 557
 558         // attempt to load index data
 559         HasIndex = LoadIndexData();
 560         success &= HasIndex;
 561     }
 562
 563     // return success status
 564     return success;
 565 }
 566
 567 bool Fasta::FastaPrivate::Rewind(void) {
 568     if ( !IsOpen ) return false;
 569     return ( fseeko(Stream, 0, SEEK_SET) == 0 );
 570 }
 571
 572 bool Fasta::FastaPrivate::WriteIndexData(void) {
 573
 574     // skip if no index file available
 575     if ( !IsIndexOpen ) return false;
 576
 577     // iterate over index entries
 578     bool success = true;
 579     stringstream indexBuffer;
 580     vector<FastaIndexData>::const_iterator indexIter = Index.begin();
 581     vector<FastaIndexData>::const_iterator indexEnd  = Index.end();
 582     for ( ; indexIter != indexEnd; ++indexIter ) {
 583
 584         // clear stream
 585         indexBuffer.str("");
 586
 587         // write data to stream
 588         const FastaIndexData& data = (*indexIter);
 589         indexBuffer << data.Name << "\t"
 590                     << data.Length << "\t"
 591                     << data.Offset << "\t"
 592                     << data.LineLength << "\t"
 593                     << data.ByteLength << endl;
 594
 595         // write stream to file
 596         success &= ( fputs(indexBuffer.str().c_str(), IndexStream) >= 0 );
 597     }
 598
 599     // return success status
 600     return success;
 601 }
 602
 603 // --------------------------------
 604 // Fasta implementation
 605
 606 Fasta::Fasta(void) {
 607     d = new FastaPrivate;
 608 }
 609
 610 Fasta::~Fasta(void) {
 611     delete d;
 612     d = 0;
 613 }
 614
 615 bool Fasta::Close(void) {
 616     return d->Close();
 617 }
 618
 619 bool Fasta::CreateIndex(const string& indexFilename) {
 620     return d->CreateIndex(indexFilename);
 621 }
 622
 623 bool Fasta::GetBase(const int& refId, const int& position, char& base) {
 624     return d->GetBase(refId, position, base);
 625 }
 626
 627 bool Fasta::GetSequence(const int& refId, const int& start, const int& stop, string& sequence) {
 628     return d->GetSequence(refId, start, stop, sequence);
 629 }
 630
 631 bool Fasta::Open(const string& filename, const string& indexFilename) {
 632     return d->Open(filename, indexFilename);
 633 }