bamtools_fasta.cpp

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