src/api/internal/BamToolsIndex_p.cpp

   1 // ***************************************************************************
   2 // BamToolsIndex.cpp (c) 2010 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 5 April 2011 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides index operations for the BamTools index format (".bti")
   9 // ***************************************************************************
  10
  11 #include <api/BamAlignment.h>
  12 #include <api/internal/BamReader_p.h>
  13 #include <api/internal/BamToolsIndex_p.h>
  14 #include <api/internal/BgzfStream_p.h>
  15 using namespace BamTools;
  16 using namespace BamTools::Internal;
  17
  18 #include <cstdio>
  19 #include <cstdlib>
  20 #include <cstring>
  21 #include <algorithm>
  22 #include <iostream>
  23 #include <map>
  24 using namespace std;
  25
  26 // static BamToolsIndex constants
  27 const int BamToolsIndex::DEFAULT_BLOCK_LENGTH = 1000;
  28 const string BamToolsIndex::BTI_EXTENSION     = ".bti";
  29 const char* const BamToolsIndex::BTI_MAGIC    = "BTI\1";
  30 const int BamToolsIndex::SIZEOF_BLOCK         = sizeof(int32_t)*2 + sizeof(int64_t);
  31
  32 // ctor
  33 BamToolsIndex::BamToolsIndex(Internal::BamReaderPrivate* reader)
  34     : BamIndex(reader)
  35     , m_indexStream(0)
  36     , m_cacheMode(BamIndex::LimitedIndexCaching)
  37     , m_blockSize(BamToolsIndex::DEFAULT_BLOCK_LENGTH)
  38     , m_inputVersion(0)
  39     , m_outputVersion(BTI_1_2) // latest version - used for writing new index files
  40 {
  41     m_isBigEndian = BamTools::SystemIsBigEndian();
  42 }
  43
  44 // dtor
  45 BamToolsIndex::~BamToolsIndex(void) {
  46     CloseFile();
  47 }
  48
  49 bool BamToolsIndex::CheckMagicNumber(void) {
  50
  51     // check 'magic number' to see if file is BTI index
  52     char magic[4];
  53     size_t elementsRead = fread(magic, sizeof(char), 4, m_indexStream);
  54     if ( elementsRead != 4 ) {
  55         cerr << "BamToolsIndex ERROR: could not read format 'magic' number" << endl;
  56         return false;
  57     }
  58
  59     if ( strncmp(magic, BamToolsIndex::BTI_MAGIC, 4) != 0 ) {
  60         cerr << "BamToolsIndex ERROR: invalid format" << endl;
  61         return false;
  62     }
  63
  64     // otherwise ok
  65     return true;
  66 }
  67
  68 // check index file version, return true if OK
  69 bool BamToolsIndex::CheckVersion(void) {
  70
  71     // read version from file
  72     size_t elementsRead = fread(&m_inputVersion, sizeof(m_inputVersion), 1, m_indexStream);
  73     if ( elementsRead != 1 ) return false;
  74     if ( m_isBigEndian ) SwapEndian_32(m_inputVersion);
  75
  76     // if version is negative, or zero
  77     if ( m_inputVersion <= 0 ) {
  78         cerr << "BamToolsIndex ERROR: could not load index file: invalid version."
  79              << endl;
  80         return false;
  81     }
  82
  83     // if version is newer than can be supported by this version of bamtools
  84     else if ( m_inputVersion > m_outputVersion ) {
  85         cerr << "BamToolsIndex ERROR: could not load index file. This version of BamTools does not recognize new index file version"
  86              << endl
  87              << "Please update BamTools to a more recent version to support this index file."
  88              << endl;
  89         return false;
  90     }
  91
  92     // ------------------------------------------------------------------
  93     // check for deprecated, unsupported versions
  94     // (typically whose format did not accomodate a particular bug fix)
  95
  96     else if ( (Version)m_inputVersion == BamToolsIndex::BTI_1_0 ) {
  97         cerr << "BamToolsIndex ERROR: could not load index file. This version of the index contains a bug related to accessing data near reference ends."
  98              << endl << endl
  99              << "Please run 'bamtools index -bti -in yourData.bam' to generate an up-to-date, fixed BTI file."
 100              << endl << endl;
 101         return false;
 102     }
 103
 104     else if ( (Version)m_inputVersion == BamToolsIndex::BTI_1_1 ) {
 105         cerr << "BamToolsIndex ERROR: could not load index file. This version of the index contains a bug related to handling empty references."
 106              << endl << endl
 107              << "Please run 'bamtools index -bti -in yourData.bam' to generate an up-to-date, fixed BTI file."
 108              << endl << endl;
 109         return false;
 110     }
 111
 112     // otherwise ok
 113     else return true;
 114 }
 115
 116 void BamToolsIndex::ClearReferenceEntry(BtiReferenceEntry& refEntry) {
 117     refEntry.ID = -1;
 118     refEntry.Blocks.clear();
 119 }
 120
 121 void BamToolsIndex::CloseFile(void) {
 122     if ( IsFileOpen() )
 123         fclose(m_indexStream);
 124     m_indexFileSummary.clear();
 125 }
 126
 127 // builds index from associated BAM file & writes out to index file
 128 bool BamToolsIndex::Create(void) {
 129
 130     // return false if BamReader is invalid or not open
 131     if ( m_reader == 0 || !m_reader->IsOpen() ) {
 132         cerr << "BamToolsIndex ERROR: BamReader is not open"
 133              << ", aborting index creation" << endl;
 134         return false;
 135     }
 136
 137     // rewind BamReader
 138     if ( !m_reader->Rewind() ) {
 139         cerr << "BamToolsIndex ERROR: could not rewind BamReader to create index"
 140              << ", aborting index creation" << endl;
 141         return false;
 142     }
 143
 144     // open new index file (read & write)
 145     string indexFilename = m_reader->Filename() + Extension();
 146     if ( !OpenFile(indexFilename, "w+b") ) {
 147         cerr << "BamStandardIndex ERROR: could not open ouput index file " << indexFilename
 148              << ", aborting index creation" << endl;
 149         return false;
 150     }
 151
 152     // initialize BtiFileSummary with number of references
 153     const int& numReferences = m_reader->GetReferenceCount();
 154     InitializeFileSummary(numReferences);
 155
 156     // initialize output file
 157     bool createdOk = true;
 158     createdOk &= WriteHeader();
 159
 160     // index building markers
 161     int32_t currentBlockCount      = 0;
 162     int64_t currentAlignmentOffset = m_reader->Tell();
 163     int32_t blockRefId             = 0;
 164     int32_t blockMaxEndPosition    = 0;
 165     int64_t blockStartOffset       = currentAlignmentOffset;
 166     int32_t blockStartPosition     = -1;
 167
 168     // plow through alignments, storing index entries
 169     BamAlignment al;
 170     BtiReferenceEntry refEntry;
 171     while ( m_reader->LoadNextAlignment(al) ) {
 172
 173         // if block contains data (not the first time through) AND alignment is on a new reference
 174         if ( currentBlockCount > 0 && al.RefID != blockRefId ) {
 175
 176             // store previous BTI block data in reference entry
 177             BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 178             refEntry.Blocks.push_back(block);
 179
 180             // write reference entry, then clear
 181             createdOk &= WriteReferenceEntry(refEntry);
 182             ClearReferenceEntry(refEntry);
 183
 184             // update markers
 185             currentBlockCount   = 0;
 186             blockMaxEndPosition = al.GetEndPosition();
 187             blockStartOffset    = currentAlignmentOffset;
 188         }
 189
 190         // if beginning of block, save first alignment's refID & position
 191         if ( currentBlockCount == 0 ) {
 192             blockRefId         = al.RefID;
 193             blockStartPosition = al.Position;
 194         }
 195
 196         // increment block counter
 197         ++currentBlockCount;
 198
 199         // check end position
 200         int32_t alignmentEndPosition = al.GetEndPosition();
 201         if ( alignmentEndPosition > blockMaxEndPosition )
 202             blockMaxEndPosition = alignmentEndPosition;
 203
 204         // if block is full, get offset for next block, reset currentBlockCount
 205         if ( currentBlockCount == m_blockSize ) {
 206
 207             // store previous block data in reference entry
 208             BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 209             refEntry.Blocks.push_back(block);
 210
 211             // update markers
 212             blockStartOffset  = m_reader->Tell();
 213             currentBlockCount = 0;
 214         }
 215
 216         // not the best name, but for the next iteration, this value will be the offset of the *current* alignment
 217         // necessary because we won't know if this next alignment is on a new reference until we actually read it
 218         currentAlignmentOffset = m_reader->Tell();
 219     }
 220
 221     // store last BTI block data in reference entry
 222     BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 223     refEntry.Blocks.push_back(block);
 224
 225     // write last reference entry, then clear
 226     createdOk &= WriteReferenceEntry(refEntry);
 227     ClearReferenceEntry(refEntry);
 228
 229     // rewind reader & return result
 230     createdOk &= m_reader->Rewind();
 231
 232     // return result
 233     return createdOk;
 234 }
 235
 236 // returns format's file extension
 237 const std::string BamToolsIndex::Extension(void) {
 238     return BamToolsIndex::BTI_EXTENSION;
 239 }
 240
 241 bool BamToolsIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) {
 242
 243     // return false ref ID is not a valid index in file summary data
 244     if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() )
 245         return false;
 246
 247     // retrieve reference index data for left bound reference
 248     BtiReferenceEntry refEntry(region.LeftRefID);
 249     if ( !ReadReferenceEntry(refEntry) ) {
 250         cerr << "BamToolsIndex ERROR: could not retrieve index data from BTI file" << endl;
 251         return false;
 252     }
 253
 254     // TODO: can this next loop be sped up?
 255     // Binary search for overlap instead of linear search perhaps.
 256
 257     // iterate over blocks on this reference
 258     BtiBlockVector::const_iterator blockIter = refEntry.Blocks.begin();
 259     BtiBlockVector::const_iterator blockEnd  = refEntry.Blocks.end();
 260     for ( ; blockIter != blockEnd; ++blockIter ) {
 261         const BtiBlock& block = (*blockIter);
 262
 263         // store offset & break if block end overlaps region start
 264         if ( block.MaxEndPosition >= region.LeftPosition ) {
 265             offset = block.StartOffset;
 266             break;
 267         }
 268     }
 269
 270     // sets to false if blocks container is empty, or if no matching block could be found
 271     *hasAlignmentsInRegion = ( blockIter != blockEnd );
 272
 273     // return success
 274     return true;
 275 }
 276
 277 // returns whether reference has alignments or no
 278 bool BamToolsIndex::HasAlignments(const int& referenceID) const {
 279     if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() )
 280         return false;
 281     const BtiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID);
 282     return ( refSummary.NumBlocks > 0 );
 283 }
 284
 285 void BamToolsIndex::InitializeFileSummary(const int& numReferences) {
 286     m_indexFileSummary.clear();
 287     for ( int i = 0; i < numReferences; ++i )
 288         m_indexFileSummary.push_back( BtiReferenceSummary() );
 289 }
 290
 291 bool BamToolsIndex::IsFileOpen(void) const {
 292     return ( m_indexStream != 0 );
 293 }
 294
 295 // attempts to use index data to jump to @region, returns success/fail
 296 // a "successful" jump indicates no error, but not whether this region has data
 297 //   * thus, the method sets a flag to indicate whether there are alignments
 298 //     available after the jump position
 299 bool BamToolsIndex::Jump(const BamTools::BamRegion& region, bool* hasAlignmentsInRegion) {
 300
 301     // clear flag
 302     *hasAlignmentsInRegion = false;
 303
 304     // skip if invalid reader or not open
 305     if ( m_reader == 0 || !m_reader->IsOpen() )
 306         return false;
 307
 308     // make sure left-bound position is valid
 309     const RefVector& references = m_reader->GetReferenceData();
 310     if ( region.LeftPosition > references.at(region.LeftRefID).RefLength )
 311         return false;
 312
 313     // calculate nearest offset to jump to
 314     int64_t offset;
 315     if ( !GetOffset(region, offset, hasAlignmentsInRegion) ) {
 316         cerr << "BamToolsIndex ERROR: could not jump"
 317              << ", unable to calculate offset for specified region" << endl;
 318         return false;
 319     }
 320
 321     // return success/failure of seek
 322     return m_reader->Seek(offset);
 323 }
 324
 325 // loads existing data from file into memory
 326 bool BamToolsIndex::Load(const std::string& filename) {
 327
 328     // attempt open index file (read-only)
 329     if ( !OpenFile(filename, "rb") ) {
 330         cerr << "BamStandardIndex ERROR: could not open input index file " << filename
 331              << ", aborting index load" << endl;
 332         return false;
 333     }
 334
 335     // attempt to load & validate BTI header data
 336     if ( !LoadHeader() ) {
 337         CloseFile();
 338         return false;
 339     }
 340
 341     // attempt to load index file summary, return success/failure
 342     return LoadFileSummary();
 343 }
 344
 345 bool BamToolsIndex::LoadFileSummary(void) {
 346
 347     // load number of reference sequences
 348     int numReferences;
 349     if ( !LoadNumReferences(numReferences) )
 350         return false;
 351
 352     // initialize file summary data
 353     InitializeFileSummary(numReferences);
 354
 355     // iterate over reference entries
 356     bool loadedOk = true;
 357     BtiFileSummary::iterator summaryIter = m_indexFileSummary.begin();
 358     BtiFileSummary::iterator summaryEnd  = m_indexFileSummary.end();
 359     for ( ; summaryIter != summaryEnd; ++summaryIter )
 360         loadedOk &= LoadReferenceSummary(*summaryIter);
 361
 362     // return result
 363     return loadedOk;
 364 }
 365
 366 bool BamToolsIndex::LoadHeader(void) {
 367
 368     // if invalid format 'magic number'
 369     if ( !CheckMagicNumber() )
 370         return false;
 371
 372     // if invalid BTI version
 373     if ( !CheckVersion() )
 374         return false;
 375
 376     // use file's BTI block size to set member variable
 377     size_t elementsRead = fread(&m_blockSize, sizeof(m_blockSize), 1, m_indexStream);
 378     if ( m_isBigEndian ) SwapEndian_32(m_blockSize);
 379     return ( elementsRead == 1 );
 380 }
 381
 382 bool BamToolsIndex::LoadNumBlocks(int& numBlocks) {
 383     size_t elementsRead = 0;
 384     elementsRead += fread(&numBlocks, sizeof(numBlocks), 1, m_indexStream);
 385     if ( m_isBigEndian ) SwapEndian_32(numBlocks);
 386     return ( elementsRead == 1 );
 387 }
 388
 389 bool BamToolsIndex::LoadNumReferences(int& numReferences) {
 390     size_t elementsRead = 0;
 391     elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
 392     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 393     return ( elementsRead == 1 );
 394 }
 395
 396 bool BamToolsIndex::LoadReferenceSummary(BtiReferenceSummary& refSummary) {
 397
 398     // load number of blocks
 399     int numBlocks;
 400     if ( !LoadNumBlocks(numBlocks) )
 401         return false;
 402
 403     // store block summary data for this reference
 404     refSummary.NumBlocks = numBlocks;
 405     refSummary.FirstBlockFilePosition = Tell();
 406
 407     // skip blocks in index file (and return status)
 408     return SkipBlocks(numBlocks);
 409 }
 410
 411 bool BamToolsIndex::OpenFile(const std::string& filename, const char* mode) {
 412
 413     // make sure any previous index file is closed
 414     CloseFile();
 415
 416     // attempt to open file
 417     m_indexStream = fopen(filename.c_str(), mode);
 418     return IsFileOpen();
 419 }
 420
 421 bool BamToolsIndex::ReadBlock(BtiBlock& block) {
 422
 423     // read in block data members
 424     size_t elementsRead = 0;
 425     elementsRead += fread(&block.MaxEndPosition, sizeof(block.MaxEndPosition), 1, m_indexStream);
 426     elementsRead += fread(&block.StartOffset,    sizeof(block.StartOffset),    1, m_indexStream);
 427     elementsRead += fread(&block.StartPosition,  sizeof(block.StartPosition),  1, m_indexStream);
 428
 429     // swap endian-ness if necessary
 430     if ( m_isBigEndian ) {
 431         SwapEndian_32(block.MaxEndPosition);
 432         SwapEndian_64(block.StartOffset);
 433         SwapEndian_32(block.StartPosition);
 434     }
 435
 436     // return success/failure
 437     return ( elementsRead == 3 );
 438 }
 439
 440 bool BamToolsIndex::ReadBlocks(const BtiReferenceSummary& refSummary, BtiBlockVector& blocks) {
 441
 442     // prep blocks container
 443     blocks.clear();
 444     blocks.reserve(refSummary.NumBlocks);
 445
 446     // skip to first block entry
 447     if ( !Seek( refSummary.FirstBlockFilePosition, SEEK_SET ) ) {
 448         cerr << "BamToolsIndex ERROR: could not seek to position "
 449              << refSummary.FirstBlockFilePosition << endl;
 450         return false;
 451     }
 452
 453     // read & store block entries
 454     bool readOk = true;
 455     BtiBlock block;
 456     for ( int i = 0; i < refSummary.NumBlocks; ++i ) {
 457         readOk &= ReadBlock(block);
 458         blocks.push_back(block);
 459     }
 460     return readOk;
 461 }
 462
 463 bool BamToolsIndex::ReadReferenceEntry(BtiReferenceEntry& refEntry) {
 464
 465     // return false if refId not valid index in file summary structure
 466     if ( refEntry.ID < 0 || refEntry.ID >= (int)m_indexFileSummary.size() )
 467         return false;
 468
 469     // use index summary to assist reading the reference's BTI blocks
 470     const BtiReferenceSummary& refSummary = m_indexFileSummary.at(refEntry.ID);
 471     return ReadBlocks(refSummary, refEntry.Blocks);
 472 }
 473
 474 bool BamToolsIndex::Seek(const int64_t& position, const int& origin) {
 475     return ( fseek64(m_indexStream, position, origin) == 0 );
 476 }
 477
 478 // change the index caching behavior
 479 void BamToolsIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) {
 480     m_cacheMode = mode;
 481     // do nothing else here ? cache mode will be ignored from now on, most likely
 482 }
 483
 484 bool BamToolsIndex::SkipBlocks(const int& numBlocks) {
 485     return Seek( numBlocks*BamToolsIndex::SIZEOF_BLOCK, SEEK_CUR );
 486 }
 487
 488 int64_t BamToolsIndex::Tell(void) const {
 489     return ftell64(m_indexStream);
 490 }
 491
 492 bool BamToolsIndex::WriteBlock(const BtiBlock& block) {
 493
 494     // copy entry data
 495     int32_t maxEndPosition = block.MaxEndPosition;
 496     int64_t startOffset    = block.StartOffset;
 497     int32_t startPosition  = block.StartPosition;
 498
 499     // swap endian-ness if necessary
 500     if ( m_isBigEndian ) {
 501         SwapEndian_32(maxEndPosition);
 502         SwapEndian_64(startOffset);
 503         SwapEndian_32(startPosition);
 504     }
 505
 506     // write the reference index entry
 507     size_t elementsWritten = 0;
 508     elementsWritten += fwrite(&maxEndPosition, sizeof(maxEndPosition), 1, m_indexStream);
 509     elementsWritten += fwrite(&startOffset,    sizeof(startOffset),    1, m_indexStream);
 510     elementsWritten += fwrite(&startPosition,  sizeof(startPosition),  1, m_indexStream);
 511     return ( elementsWritten == 3 );
 512 }
 513
 514 bool BamToolsIndex::WriteBlocks(const BtiBlockVector& blocks) {
 515     bool writtenOk = true;
 516     BtiBlockVector::const_iterator blockIter = blocks.begin();
 517     BtiBlockVector::const_iterator blockEnd  = blocks.end();
 518     for ( ; blockIter != blockEnd; ++blockIter )
 519         writtenOk &= WriteBlock(*blockIter);
 520     return writtenOk;
 521 }
 522
 523 bool BamToolsIndex::WriteHeader(void) {
 524
 525     size_t elementsWritten = 0;
 526
 527     // write BTI index format 'magic number'
 528     elementsWritten += fwrite(BamToolsIndex::BTI_MAGIC, 1, 4, m_indexStream);
 529
 530     // write BTI index format version
 531     int32_t currentVersion = (int32_t)m_outputVersion;
 532     if ( m_isBigEndian ) SwapEndian_32(currentVersion);
 533     elementsWritten += fwrite(&currentVersion, sizeof(currentVersion), 1, m_indexStream);
 534
 535     // write block size
 536     int32_t blockSize = m_blockSize;
 537     if ( m_isBigEndian ) SwapEndian_32(blockSize);
 538     elementsWritten += fwrite(&blockSize, sizeof(blockSize), 1, m_indexStream);
 539
 540     // write number of references
 541     int32_t numReferences = m_indexFileSummary.size();
 542     cerr << "Writing numReferences = " << numReferences << endl;
 543     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 544     elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, m_indexStream);
 545
 546     // return success/failure of write
 547     return ( elementsWritten == 7 );
 548 }
 549
 550 bool BamToolsIndex::WriteReferenceEntry(const BtiReferenceEntry& refEntry) {
 551
 552     size_t elementsWritten = 0;
 553
 554     // write number of blocks this reference
 555     uint32_t numBlocks = refEntry.Blocks.size();
 556     if ( m_isBigEndian ) SwapEndian_32(numBlocks);
 557     elementsWritten += fwrite(&numBlocks, sizeof(numBlocks), 1, m_indexStream);
 558
 559     // write actual block entries
 560     const bool blocksOk = WriteBlocks(refEntry.Blocks);
 561
 562     // return success/fail
 563     return ( elementsWritten == 1) && blocksOk;
 564 }