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