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 <cstdio>
  12 #include <cstdlib>
  13 #include <cstring>
  14 #include <fstream>
  15 #include <iostream>
  16 #include <sstream>
  17 #include <vector>
  18 #include "bamtools_fasta.h"
  19 using namespace std;
  20 using namespace BamTools;
  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         // seek to beginning of sequence data
 250         if ( fseeko(Stream, referenceData.Offset, SEEK_SET) != 0 ) {
 251             cerr << "FASTA error : could not sek in file" << endl;
 252             return false;
 253         }
 254
 255         // retrieve sequence
 256         string sequence = "";
 257         if ( !GetNextSequence(sequence) ) {
 258             cerr << "FASTA error : could not retrieve base from FASTA file" << endl;
 259             return false;
 260         }
 261
 262         // set base & return success
 263         base = sequence.at(position);
 264         return true;
 265     }
 266
 267     // else plow through sequentially
 268     else {
 269
 270         // rewind FASTA file
 271         if ( !Rewind() ) {
 272             cerr << "FASTA error : could not rewind FASTA file" << endl;
 273             return false;
 274         }
 275
 276         // iterate through fasta entries
 277         int currentId = 0;
 278         string header = "";
 279         string sequence = "";
 280
 281         // get first entry
 282         GetNextHeader(header);
 283         GetNextSequence(sequence);
 284
 285         while ( currentId != refId ) {
 286             GetNextHeader(header);
 287             GetNextSequence(sequence);
 288             ++currentId;
 289         }
 290
 291         // get desired base from sequence
 292         // TODO: error reporting on invalid position
 293         if ( currentId == refId && (sequence.length() >= (size_t)position) ) {
 294             base = sequence.at(position);
 295             return true;
 296         }
 297
 298         // could not get sequence
 299         return false;
 300     }
 301
 302     // return success
 303     return true;
 304 }
 305
 306 bool Fasta::FastaPrivate::GetNameFromHeader(const string& header, string& name) {
 307
 308     // get rid of the leading greater than sign
 309     string s = header.substr(1);
 310
 311     // extract the first non-whitespace segment
 312     char* pName = (char*)s.data();
 313     unsigned int nameLen = (unsigned int)s.size();
 314
 315     unsigned int start = 0;
 316     while ( (pName[start] == 32) || (pName[start] == 9) || (pName[start] == 10) || (pName[start] == 13) ) {
 317         start++;
 318         if ( start == nameLen )
 319             break;
 320     }
 321
 322     unsigned int stop  = start;
 323     if ( stop < nameLen ) {
 324         while( (pName[stop] != 32) && (pName[stop] != 9) && (pName[stop] != 10) && (pName[stop] != 13) ) {
 325             stop++;
 326             if ( stop == nameLen )
 327                 break;
 328         }
 329     }
 330
 331     if ( start == stop ) {
 332         cerr << "FASTA error : could not parse read name from FASTA header" << endl;
 333         return false;
 334     }
 335
 336     name = s.substr(start, stop - start).c_str();
 337     return true;
 338 }
 339
 340 bool Fasta::FastaPrivate::GetNextHeader(string& header) {
 341
 342     // validate input stream
 343     if ( !IsOpen || feof(Stream) )
 344         return false;
 345
 346     // read in header line
 347     char buffer[1024];
 348     if ( fgets(buffer, 1024, Stream) == 0 ) {
 349         cerr << "FASTA error : could not read from file" << endl;
 350         return false;
 351     }
 352
 353     // make sure it's a FASTA header
 354     if ( buffer[0] != '>' ) {
 355         cerr << "FASTA error : expected header ('>'), instead : " << buffer[0] << endl;
 356         return false;
 357     }
 358
 359     // import buffer contents to header string
 360     stringstream headerBuffer("");
 361     headerBuffer << buffer;
 362     header = headerBuffer.str();
 363
 364     // return success
 365     return true;
 366 }
 367
 368 bool Fasta::FastaPrivate::GetNextSequence(string& sequence) {
 369
 370     // validate input stream
 371     if ( !IsOpen || feof(Stream) )
 372         return false;
 373
 374     // read in sequence
 375     char buffer[1024];
 376     ostringstream seqBuffer("");
 377     while(true) {
 378
 379         char ch = fgetc(Stream);
 380         ungetc(ch, Stream);
 381         if( (ch == '>') || feof(Stream) )
 382               break;
 383
 384         if ( fgets(buffer, 1024, Stream) == 0 ) {
 385             cerr << "FASTA error : could not read from file" << endl;
 386             return false;
 387         }
 388
 389         Chomp(buffer);
 390         seqBuffer << buffer;
 391     }
 392
 393     // import buffer contents to sequence string
 394     sequence = seqBuffer.str();
 395
 396     // return success
 397     return true;
 398 }
 399
 400 bool Fasta::FastaPrivate::GetSequence(const int& refId, const int& start, const int& stop, string& sequence) {
 401
 402     // make sure FASTA file is open
 403     if ( !IsOpen ) {
 404         cerr << "FASTA error : file not open for reading" << endl;
 405         return false;
 406     }
 407
 408     // use index if available
 409     if ( HasIndex && !Index.empty() ) {
 410
 411         // validate reference id
 412         if ( (refId < 0) || (refId >= (int)Index.size()) ) {
 413             cerr << "FASTA error: invalid refId specified: " << refId << endl;
 414             return false;
 415         }
 416
 417         // retrieve reference index data
 418         const FastaIndexData& referenceData = Index.at(refId);
 419
 420         // validate stop position
 421         if ( (start < 0) || (start > stop) || (stop > referenceData.Length) ) {
 422             cerr << "FASTA error: invalid start/stop positions specified: " << start << ", " << stop << endl;
 423             return false;
 424         }
 425
 426         // seek to beginning of sequence data
 427         if ( fseeko(Stream, referenceData.Offset, SEEK_SET) != 0 ) {
 428             cerr << "FASTA error : could not sek in file" << endl;
 429             return false;
 430         }
 431
 432         // retrieve full sequence
 433         string fullSequence = "";
 434         if ( !GetNextSequence(fullSequence) ) {
 435             cerr << "FASTA error : could not retrieve sequence from FASTA file" << endl;
 436             return false;
 437         }
 438
 439         // set sub-sequence & return success
 440         const int seqLength = (stop - start) + 1;
 441         sequence = fullSequence.substr(start, seqLength);
 442         return true;
 443     }
 444
 445     // else plow through sequentially
 446     else {
 447
 448         // rewind FASTA file
 449         if ( !Rewind() ) {
 450             cerr << "FASTA error : could not rewind FASTA file" << endl;
 451             return false;
 452         }
 453
 454         // iterate through fasta entries
 455         int currentId = 0;
 456         string header = "";
 457         string fullSequence = "";
 458
 459         // get first entry
 460         GetNextHeader(header);
 461         GetNextSequence(fullSequence);
 462
 463         while ( currentId != refId ) {
 464             GetNextHeader(header);
 465             GetNextSequence(fullSequence);
 466             ++currentId;
 467         }
 468
 469         // get desired substring from sequence
 470         // TODO: error reporting on invalid start/stop positions
 471         if ( currentId == refId && (fullSequence.length() >= (size_t)stop) ) {
 472             const int seqLength = (stop - start) + 1;
 473             sequence = fullSequence.substr(start, seqLength);
 474             return true;
 475         }
 476
 477         // could not get sequence
 478         return false;
 479     }
 480
 481     // return success
 482     return true;
 483 }
 484
 485 bool Fasta::FastaPrivate::LoadIndexData(void) {
 486
 487     // skip if no index file available
 488     if ( !IsIndexOpen ) return false;
 489
 490     // clear any prior index data
 491     Index.clear();
 492
 493     char buffer[1024];
 494     stringstream indexBuffer;
 495     while ( true ) {
 496
 497         char c = fgetc(IndexStream);
 498         if ( (c == '\n') || feof(IndexStream) ) break;
 499         ungetc(c, IndexStream);
 500
 501         // clear index buffer
 502         indexBuffer.str("");
 503
 504         // read line from index file
 505         if ( fgets(buffer, 1024, IndexStream) == 0 ) {
 506             cerr << "FASTA LoadIndexData() error : could not read from index file" << endl;
 507             HasIndex = false;
 508             return false;
 509         }
 510
 511         // store line in indexBuffer
 512         indexBuffer << buffer;
 513
 514         // retrieve fasta index data from line
 515         FastaIndexData data;
 516         indexBuffer >> data.Name;
 517         indexBuffer >> data.Length;
 518         indexBuffer >> data.Offset;
 519         indexBuffer >> data.LineLength;
 520         indexBuffer >> data.ByteLength;
 521
 522         // store index entry
 523         Index.push_back(data);
 524     }
 525
 526     return true;
 527 }
 528
 529 bool Fasta::FastaPrivate::Open(const string& filename, const string& indexFilename) {
 530
 531     bool success = true;
 532
 533     // open FASTA filename
 534     Stream = fopen(filename.c_str(), "rb");
 535     if ( !Stream ) {
 536         cerr << "FASTA error: Could not open " << filename << " for reading" << endl;
 537         return false;
 538     }
 539     IsOpen = true;
 540     success &= IsOpen;
 541
 542     // open index file if it exists
 543     if ( !indexFilename.empty() ) {
 544         IndexStream = fopen(indexFilename.c_str(), "rb");
 545         if ( !IndexStream ) {
 546             cerr << "FASTA error : Could not open " << indexFilename << " for reading." << endl;
 547             return false;
 548         }
 549         IsIndexOpen = true;
 550         success &= IsIndexOpen;
 551
 552         // attempt to load index data
 553         HasIndex = LoadIndexData();
 554         success &= HasIndex;
 555     }
 556
 557     // return success status
 558     return success;
 559 }
 560
 561 bool Fasta::FastaPrivate::Rewind(void) {
 562     if ( !IsOpen ) return false;
 563     return ( fseeko(Stream, 0, SEEK_SET) == 0 );
 564 }
 565
 566 bool Fasta::FastaPrivate::WriteIndexData(void) {
 567
 568     // skip if no index file available
 569     if ( !IsIndexOpen ) return false;
 570
 571     // iterate over index entries
 572     bool success = true;
 573     stringstream indexBuffer;
 574     vector<FastaIndexData>::const_iterator indexIter = Index.begin();
 575     vector<FastaIndexData>::const_iterator indexEnd  = Index.end();
 576     for ( ; indexIter != indexEnd; ++indexIter ) {
 577
 578         // clear stream
 579         indexBuffer.str("");
 580
 581         // write data to stream
 582         const FastaIndexData& data = (*indexIter);
 583         indexBuffer << data.Name << "\t"
 584                     << data.Length << "\t"
 585                     << data.Offset << "\t"
 586                     << data.LineLength << "\t"
 587                     << data.ByteLength << endl;
 588
 589         // write stream to file
 590         success &= ( fputs(indexBuffer.str().c_str(), IndexStream) >= 0 );
 591     }
 592
 593     // return success status
 594     return success;
 595 }
 596
 597 // --------------------------------
 598 // Fasta implementation
 599
 600 Fasta::Fasta(void) {
 601     d = new FastaPrivate;
 602 }
 603
 604 Fasta::~Fasta(void) {
 605     delete d;
 606     d = 0;
 607 }
 608
 609 bool Fasta::Close(void) {
 610     return d->Close();
 611 }
 612
 613 bool Fasta::CreateIndex(const string& indexFilename) {
 614     return d->CreateIndex(indexFilename);
 615 }
 616
 617 bool Fasta::GetBase(const int& refId, const int& position, char& base) {
 618     return d->GetBase(refId, position, base);
 619 }
 620
 621 bool Fasta::GetSequence(const int& refId, const int& start, const int& stop, string& sequence) {
 622     return d->GetSequence(refId, start, stop, sequence);
 623 }
 624
 625 bool Fasta::Open(const string& filename, const string& indexFilename) {
 626     return d->Open(filename, indexFilename);
 627 }