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