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