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: 22 November 2010 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides index operations for the BamTools index format (".bti")
   9 // ***************************************************************************
  10
  11 #include <api/BamAlignment.h>
  12 #include <api/BamReader.h>
  13 #include <api/BGZF.h>
  14 #include <api/internal/BamToolsIndex_p.h>
  15 using namespace BamTools;
  16 using namespace BamTools::Internal;
  17
  18 #include <cstdio>
  19 #include <cstdlib>
  20 #include <algorithm>
  21 #include <iostream>
  22 #include <map>
  23 using namespace std;
  24
  25 BamToolsIndex::BamToolsIndex(BgzfData* bgzf, BamReader* reader)
  26     : BamIndex(bgzf, reader)
  27     , m_blockSize(1000)
  28     , m_dataBeginOffset(0)
  29     , m_hasFullDataCache(false)
  30     , m_inputVersion(0)
  31     , m_outputVersion(BTI_1_2) // latest version - used for writing new index files
  32 {
  33     m_isBigEndian = BamTools::SystemIsBigEndian();
  34 }
  35
  36 // dtor
  37 BamToolsIndex::~BamToolsIndex(void) {
  38     ClearAllData();
  39 }
  40
  41 // creates index data (in-memory) from current reader data
  42 bool BamToolsIndex::Build(void) {
  43
  44     // be sure reader & BGZF file are valid & open for reading
  45     if ( m_reader == 0 || m_BGZF == 0 || !m_BGZF->IsOpen )
  46         return false;
  47
  48     // move file pointer to beginning of alignments
  49     if ( !m_reader->Rewind() ) return false;
  50
  51     // initialize index data structure with space for all references
  52     const int numReferences = (int)m_references.size();
  53     m_indexData.clear();
  54     m_hasFullDataCache = false;
  55     SetReferenceCount(numReferences);
  56
  57     // set up counters and markers
  58     int32_t currentBlockCount      = 0;
  59     int64_t currentAlignmentOffset = m_BGZF->Tell();
  60     int32_t blockRefId             = 0;
  61     int32_t blockMaxEndPosition    = 0;
  62     int64_t blockStartOffset       = currentAlignmentOffset;
  63     int32_t blockStartPosition     = -1;
  64
  65     // plow through alignments, storing index entries
  66     BamAlignment al;
  67     while ( m_reader->GetNextAlignmentCore(al) ) {
  68
  69         // if block contains data (not the first time through) AND alignment is on a new reference
  70         if ( currentBlockCount > 0 && al.RefID != blockRefId ) {
  71
  72             // store previous data
  73             BamToolsIndexEntry entry(blockMaxEndPosition, blockStartOffset, blockStartPosition);
  74             SaveOffsetEntry(blockRefId, entry);
  75
  76             // intialize new block for current alignment's reference
  77             currentBlockCount   = 0;
  78             blockMaxEndPosition = al.GetEndPosition();
  79             blockStartOffset    = currentAlignmentOffset;
  80         }
  81
  82         // if beginning of block, save first alignment's refID & position
  83         if ( currentBlockCount == 0 ) {
  84             blockRefId         = al.RefID;
  85             blockStartPosition = al.Position;
  86         }
  87
  88         // increment block counter
  89         ++currentBlockCount;
  90
  91         // check end position
  92         int32_t alignmentEndPosition = al.GetEndPosition();
  93         if ( alignmentEndPosition > blockMaxEndPosition )
  94             blockMaxEndPosition = alignmentEndPosition;
  95
  96         // if block is full, get offset for next block, reset currentBlockCount
  97         if ( currentBlockCount == m_blockSize ) {
  98             BamToolsIndexEntry entry(blockMaxEndPosition, blockStartOffset, blockStartPosition);
  99             SaveOffsetEntry(blockRefId, entry);
 100             blockStartOffset  = m_BGZF->Tell();
 101             currentBlockCount = 0;
 102         }
 103
 104         // not the best name, but for the next iteration, this value will be the offset of the *current* alignment
 105         // necessary because we won't know if this next alignment is on a new reference until we actually read it
 106         currentAlignmentOffset = m_BGZF->Tell();
 107     }
 108
 109     // store final block with data
 110     BamToolsIndexEntry entry(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 111     SaveOffsetEntry(blockRefId, entry);
 112
 113     // set flag
 114     m_hasFullDataCache = true;
 115
 116     // return success/failure of rewind
 117     return m_reader->Rewind();
 118 }
 119
 120 // check index file magic number, return true if OK
 121 bool BamToolsIndex::CheckMagicNumber(void) {
 122
 123     // see if index is valid BAM index
 124     char magic[4];
 125     size_t elementsRead = fread(magic, 1, 4, m_indexStream);
 126     if ( elementsRead != 4 ) return false;
 127     if ( strncmp(magic, "BTI\1", 4) != 0 ) {
 128         fprintf(stderr, "Problem with index file - invalid format.\n");
 129         return false;
 130     }
 131
 132     // otherwise ok
 133     return true;
 134 }
 135
 136 // check index file version, return true if OK
 137 bool BamToolsIndex::CheckVersion(void) {
 138
 139     // read version from file
 140     size_t elementsRead = fread(&m_inputVersion, sizeof(m_inputVersion), 1, m_indexStream);
 141     if ( elementsRead != 1 ) return false;
 142     if ( m_isBigEndian ) SwapEndian_32(m_inputVersion);
 143
 144     // if version is negative, or zero
 145     if ( m_inputVersion <= 0 ) {
 146         fprintf(stderr, "Problem with index file - invalid version.\n");
 147         return false;
 148     }
 149
 150     // if version is newer than can be supported by this version of bamtools
 151     else if ( m_inputVersion > m_outputVersion ) {
 152         fprintf(stderr, "Problem with index file - attempting to use an outdated version of BamTools with a newer index file.\n");
 153         fprintf(stderr, "Please update BamTools to a more recent version to support this index file.\n");
 154         return false;
 155     }
 156
 157     // ------------------------------------------------------------------
 158     // check for deprecated, unsupported versions
 159     // (typically whose format did not accomodate a particular bug fix)
 160
 161     else if ( (Version)m_inputVersion == BTI_1_0 ) {
 162         fprintf(stderr, "\nProblem with index file - this version of the index contains a bug related to accessing data near reference ends.\n");
 163         fprintf(stderr, "\nPlease run \'bamtools index -bti -in yourData.bam\' to generate an up-to-date BamToolsIndex.\n\n");
 164         return false;
 165     }
 166
 167     else if ( (Version)m_inputVersion == BTI_1_1 ) {
 168         fprintf(stderr, "\nProblem with index file - this version of the index contains a bug related to handling empty references.\n");
 169         fprintf(stderr, "\nPlease run \'bamtools index -bti -in yourData.bam\' to generate an up-to-date BamToolsIndex.\n\n");
 170         return false;
 171     }
 172
 173     // otherwise ok
 174     else return true;
 175 }
 176
 177 // clear all current index offset data in memory
 178 void BamToolsIndex::ClearAllData(void) {
 179     BamToolsIndexData::const_iterator indexIter = m_indexData.begin();
 180     BamToolsIndexData::const_iterator indexEnd  = m_indexData.end();
 181     for ( ; indexIter != indexEnd; ++indexIter ) {
 182         const int& refId = (*indexIter).first;
 183         ClearReferenceOffsets(refId);
 184     }
 185 }
 186
 187 // clear all index offset data for desired reference
 188 void BamToolsIndex::ClearReferenceOffsets(const int& refId) {
 189     if ( m_indexData.find(refId) == m_indexData.end() ) return;
 190     vector<BamToolsIndexEntry>& offsets = m_indexData[refId].Offsets;
 191     offsets.clear();
 192     m_hasFullDataCache = false;
 193 }
 194
 195 // return file position after header metadata
 196 const off_t BamToolsIndex::DataBeginOffset(void) const {
 197     return m_dataBeginOffset;
 198 }
 199
 200 // calculate BAM file offset for desired region
 201 // return true if no error (*NOT* equivalent to "has alignments or valid offset")
 202 //   check @hasAlignmentsInRegion to determine this status
 203 // @region - target region
 204 // @offset - resulting seek target
 205 // @hasAlignmentsInRegion - sometimes a file just lacks data in region, this flag indicates that status
 206 // N.B. - ignores isRightBoundSpecified
 207 bool BamToolsIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) {
 208
 209     // return false if leftBound refID is not found in index data
 210     BamToolsIndexData::const_iterator indexIter = m_indexData.find(region.LeftRefID);
 211     if ( indexIter == m_indexData.end()) return false;
 212
 213     // load index data for region if not already cached
 214     if ( !IsDataLoaded(region.LeftRefID) ) {
 215         bool loadedOk = true;
 216         loadedOk &= SkipToReference(region.LeftRefID);
 217         loadedOk &= LoadReference(region.LeftRefID);
 218         if ( !loadedOk ) return false;
 219     }
 220
 221     // localize index data for this reference (& sanity check that data actually exists)
 222     indexIter = m_indexData.find(region.LeftRefID);
 223     if ( indexIter == m_indexData.end()) return false;
 224     const vector<BamToolsIndexEntry>& referenceOffsets = (*indexIter).second.Offsets;
 225     if ( referenceOffsets.empty() ) return false;
 226
 227     // -------------------------------------------------------
 228     // calculate nearest index to jump to
 229
 230     // save first offset
 231     offset = (*referenceOffsets.begin()).StartOffset;
 232
 233     // iterate over offsets entries on this reference
 234     vector<BamToolsIndexEntry>::const_iterator offsetIter = referenceOffsets.begin();
 235     vector<BamToolsIndexEntry>::const_iterator offsetEnd  = referenceOffsets.end();
 236     for ( ; offsetIter != offsetEnd; ++offsetIter ) {
 237         const BamToolsIndexEntry& entry = (*offsetIter);
 238         // break if alignment 'entry' overlaps region
 239         if ( entry.MaxEndPosition >= region.LeftPosition ) break;
 240         offset = (*offsetIter).StartOffset;
 241     }
 242
 243     // set flag based on whether an index entry was found for this region
 244     *hasAlignmentsInRegion = ( offsetIter != offsetEnd );
 245
 246     // if cache mode set to none, dump the data we just loaded
 247     if (m_cacheMode == BamIndex::NoIndexCaching )
 248         ClearReferenceOffsets(region.LeftRefID);
 249
 250     // return success
 251     return true;
 252 }
 253
 254 // returns whether reference has alignments or no
 255 bool BamToolsIndex::HasAlignments(const int& refId) const {
 256
 257     BamToolsIndexData::const_iterator indexIter = m_indexData.find(refId);
 258     if ( indexIter == m_indexData.end()) return false;
 259     const BamToolsReferenceEntry& refEntry = (*indexIter).second;
 260     return refEntry.HasAlignments;
 261 }
 262
 263 // return true if all index data is cached
 264 bool BamToolsIndex::HasFullDataCache(void) const {
 265     return m_hasFullDataCache;
 266 }
 267
 268 // returns true if index cache has data for desired reference
 269 bool BamToolsIndex::IsDataLoaded(const int& refId) const {
 270
 271     BamToolsIndexData::const_iterator indexIter = m_indexData.find(refId);
 272     if ( indexIter == m_indexData.end()) return false;
 273     const BamToolsReferenceEntry& refEntry = (*indexIter).second;
 274
 275     if ( !refEntry.HasAlignments ) return true; // no data period
 276
 277     // return whether offsets list contains data
 278     return !refEntry.Offsets.empty();
 279 }
 280
 281 // attempts to use index to jump to region; returns success/fail
 282 bool BamToolsIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
 283
 284     // clear flag
 285     *hasAlignmentsInRegion = false;
 286
 287     // check valid BamReader state
 288     if ( m_reader == 0 || m_BGZF == 0 || !m_reader->IsOpen() ) {
 289         fprintf(stderr, "ERROR: Could not jump: invalid BamReader state.\n");
 290         return false;
 291     }
 292
 293     // make sure left-bound position is valid
 294     if ( region.LeftPosition > m_references.at(region.LeftRefID).RefLength )
 295         return false;
 296
 297     // calculate nearest offset to jump to
 298     int64_t offset;
 299     if ( !GetOffset(region, offset, hasAlignmentsInRegion) ) {
 300         fprintf(stderr, "ERROR: Could not jump - unable to calculate offset for specified region.\n");
 301         return false;
 302     }
 303
 304     // return success/failure of seek
 305     return m_BGZF->Seek(offset);
 306 }
 307
 308 // clears index data from all references except the first
 309 void BamToolsIndex::KeepOnlyFirstReferenceOffsets(void) {
 310     BamToolsIndexData::const_iterator indexBegin = m_indexData.begin();
 311     KeepOnlyReferenceOffsets( (*indexBegin).first );
 312 }
 313
 314 // clears index data from all references except the one specified
 315 void BamToolsIndex::KeepOnlyReferenceOffsets(const int& refId) {
 316     BamToolsIndexData::iterator mapIter = m_indexData.begin();
 317     BamToolsIndexData::iterator mapEnd  = m_indexData.end();
 318     for ( ; mapIter != mapEnd; ++mapIter ) {
 319         const int entryRefId = (*mapIter).first;
 320         if ( entryRefId != refId )
 321             ClearReferenceOffsets(entryRefId);
 322     }
 323 }
 324
 325 // load index data for all references, return true if loaded OK
 326 bool BamToolsIndex::LoadAllReferences(bool saveData) {
 327
 328     // skip if data already loaded
 329     if ( m_hasFullDataCache ) return true;
 330
 331     // read in number of references
 332     int32_t numReferences;
 333     if ( !LoadReferenceCount(numReferences) ) return false;
 334     //SetReferenceCount(numReferences);
 335
 336     // iterate over reference entries
 337     bool loadedOk = true;
 338     for ( int i = 0; i < numReferences; ++i )
 339         loadedOk &= LoadReference(i, saveData);
 340
 341     // set flag
 342     if ( loadedOk && saveData )
 343         m_hasFullDataCache = true;
 344
 345     // return success/failure of load
 346     return loadedOk;
 347 }
 348
 349 // load header data from index file, return true if loaded OK
 350 bool BamToolsIndex::LoadHeader(void) {
 351
 352     // check magic number
 353     if ( !CheckMagicNumber() ) return false;
 354
 355     // check BTI version
 356     if ( !CheckVersion() ) return false;
 357
 358     // read in block size
 359     size_t elementsRead = fread(&m_blockSize, sizeof(m_blockSize), 1, m_indexStream);
 360     if ( elementsRead != 1 ) return false;
 361     if ( m_isBigEndian ) SwapEndian_32(m_blockSize);
 362
 363     // store offset of beginning of data
 364     m_dataBeginOffset = ftell64(m_indexStream);
 365
 366     // return success/failure of load
 367     return (elementsRead == 1);
 368 }
 369
 370 // load a single index entry from file, return true if loaded OK
 371 // @saveData - save data in memory if true, just read & discard if false
 372 bool BamToolsIndex::LoadIndexEntry(const int& refId, bool saveData) {
 373
 374     // read in index entry data members
 375     size_t elementsRead = 0;
 376     BamToolsIndexEntry entry;
 377     elementsRead += fread(&entry.MaxEndPosition, sizeof(entry.MaxEndPosition), 1, m_indexStream);
 378     elementsRead += fread(&entry.StartOffset,    sizeof(entry.StartOffset),    1, m_indexStream);
 379     elementsRead += fread(&entry.StartPosition,  sizeof(entry.StartPosition),  1, m_indexStream);
 380     if ( elementsRead != 3 ) {
 381         cerr << "Error reading index entry. Expected 3 elements, read in: " << elementsRead << endl;
 382         return false;
 383     }
 384
 385     // swap endian-ness if necessary
 386     if ( m_isBigEndian ) {
 387         SwapEndian_32(entry.MaxEndPosition);
 388         SwapEndian_64(entry.StartOffset);
 389         SwapEndian_32(entry.StartPosition);
 390     }
 391
 392     // save data
 393     if ( saveData )
 394         SaveOffsetEntry(refId, entry);
 395
 396     // return success/failure of load
 397     return true;
 398 }
 399
 400 // load a single reference from file, return true if loaded OK
 401 // @saveData - save data in memory if true, just read & discard if false
 402 bool BamToolsIndex::LoadFirstReference(bool saveData) {
 403     BamToolsIndexData::const_iterator indexBegin = m_indexData.begin();
 404     return LoadReference( (*indexBegin).first, saveData );
 405 }
 406
 407 // load a single reference from file, return true if loaded OK
 408 // @saveData - save data in memory if true, just read & discard if false
 409 bool BamToolsIndex::LoadReference(const int& refId, bool saveData) {
 410
 411     // read in number of offsets for this reference
 412     uint32_t numOffsets;
 413     size_t elementsRead = fread(&numOffsets, sizeof(numOffsets), 1, m_indexStream);
 414     if ( elementsRead != 1 ) return false;
 415     if ( m_isBigEndian ) SwapEndian_32(numOffsets);
 416
 417     // initialize offsets container for this reference
 418     SetOffsetCount(refId, (int)numOffsets);
 419
 420     // iterate over offset entries
 421     for ( unsigned int j = 0; j < numOffsets; ++j )
 422         LoadIndexEntry(refId, saveData);
 423
 424     // return success/failure of load
 425     return true;
 426 }
 427
 428 // loads number of references, return true if loaded OK
 429 bool BamToolsIndex::LoadReferenceCount(int& numReferences) {
 430
 431     size_t elementsRead = 0;
 432
 433     // read reference count
 434     elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
 435     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 436
 437     // return success/failure of load
 438     return ( elementsRead == 1 );
 439 }
 440
 441 // saves an index offset entry in memory
 442 void BamToolsIndex::SaveOffsetEntry(const int& refId, const BamToolsIndexEntry& entry) {
 443     BamToolsReferenceEntry& refEntry = m_indexData[refId];
 444     refEntry.HasAlignments = true;
 445     refEntry.Offsets.push_back(entry);
 446 }
 447
 448 // pre-allocates size for offset vector
 449 void BamToolsIndex::SetOffsetCount(const int& refId, const int& offsetCount) {
 450     BamToolsReferenceEntry& refEntry = m_indexData[refId];
 451     refEntry.Offsets.reserve(offsetCount);
 452     refEntry.HasAlignments = ( offsetCount > 0);
 453 }
 454
 455 // initializes index data structure to hold @count references
 456 void BamToolsIndex::SetReferenceCount(const int& count) {
 457     for ( int i = 0; i < count; ++i )
 458         m_indexData[i].HasAlignments = false;
 459 }
 460
 461 // position file pointer to first reference begin, return true if skipped OK
 462 bool BamToolsIndex::SkipToFirstReference(void) {
 463     BamToolsIndexData::const_iterator indexBegin = m_indexData.begin();
 464     return SkipToReference( (*indexBegin).first );
 465 }
 466
 467 // position file pointer to desired reference begin, return true if skipped OK
 468 bool BamToolsIndex::SkipToReference(const int& refId) {
 469
 470     // attempt rewind
 471     if ( !Rewind() ) return false;
 472
 473     // read in number of references
 474     int32_t numReferences;
 475     size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
 476     if ( elementsRead != 1 ) return false;
 477     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 478
 479     // iterate over reference entries
 480     bool skippedOk = true;
 481     int currentRefId = 0;
 482     while (currentRefId != refId) {
 483         skippedOk &= LoadReference(currentRefId, false);
 484         ++currentRefId;
 485     }
 486
 487     // return success/failure of skip
 488     return skippedOk;
 489 }
 490
 491 // write header to new index file
 492 bool BamToolsIndex::WriteHeader(void) {
 493
 494     size_t elementsWritten = 0;
 495
 496     // write BTI index format 'magic number'
 497     elementsWritten += fwrite("BTI\1", 1, 4, m_indexStream);
 498
 499     // write BTI index format version
 500     int32_t currentVersion = (int32_t)m_outputVersion;
 501     if ( m_isBigEndian ) SwapEndian_32(currentVersion);
 502     elementsWritten += fwrite(&currentVersion, sizeof(currentVersion), 1, m_indexStream);
 503
 504     // write block size
 505     int32_t blockSize = m_blockSize;
 506     if ( m_isBigEndian ) SwapEndian_32(blockSize);
 507     elementsWritten += fwrite(&blockSize, sizeof(blockSize), 1, m_indexStream);
 508
 509     // store offset of beginning of data
 510     m_dataBeginOffset = ftell64(m_indexStream);
 511
 512     // return success/failure of write
 513     return ( elementsWritten == 6 );
 514 }
 515
 516 // write index data for all references to new index file
 517 bool BamToolsIndex::WriteAllReferences(void) {
 518
 519     size_t elementsWritten = 0;
 520
 521     // write number of references
 522     int32_t numReferences = (int32_t)m_indexData.size();
 523     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 524     elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, m_indexStream);
 525
 526     // iterate through references in index
 527     bool refOk = true;
 528     BamToolsIndexData::const_iterator refIter = m_indexData.begin();
 529     BamToolsIndexData::const_iterator refEnd  = m_indexData.end();
 530     for ( ; refIter != refEnd; ++refIter )
 531         refOk &= WriteReferenceEntry( (*refIter).second );
 532
 533     return ( (elementsWritten == 1) && refOk );
 534 }
 535
 536 // write current reference index data to new index file
 537 bool BamToolsIndex::WriteReferenceEntry(const BamToolsReferenceEntry& refEntry) {
 538
 539     size_t elementsWritten = 0;
 540
 541     // write number of offsets listed for this reference
 542     uint32_t numOffsets = refEntry.Offsets.size();
 543     if ( m_isBigEndian ) SwapEndian_32(numOffsets);
 544     elementsWritten += fwrite(&numOffsets, sizeof(numOffsets), 1, m_indexStream);
 545
 546     // iterate over offset entries
 547     bool entriesOk = true;
 548     vector<BamToolsIndexEntry>::const_iterator offsetIter = refEntry.Offsets.begin();
 549     vector<BamToolsIndexEntry>::const_iterator offsetEnd  = refEntry.Offsets.end();
 550     for ( ; offsetIter != offsetEnd; ++offsetIter )
 551         entriesOk &= WriteIndexEntry( (*offsetIter) );
 552
 553     return ( (elementsWritten == 1) && entriesOk );
 554 }
 555
 556 // write current index offset entry to new index file
 557 bool BamToolsIndex::WriteIndexEntry(const BamToolsIndexEntry& entry) {
 558
 559     // copy entry data
 560     int32_t maxEndPosition = entry.MaxEndPosition;
 561     int64_t startOffset    = entry.StartOffset;
 562     int32_t startPosition  = entry.StartPosition;
 563
 564     // swap endian-ness if necessary
 565     if ( m_isBigEndian ) {
 566         SwapEndian_32(maxEndPosition);
 567         SwapEndian_64(startOffset);
 568         SwapEndian_32(startPosition);
 569     }
 570
 571     // write the reference index entry
 572     size_t elementsWritten = 0;
 573     elementsWritten += fwrite(&maxEndPosition, sizeof(maxEndPosition), 1, m_indexStream);
 574     elementsWritten += fwrite(&startOffset,    sizeof(startOffset),    1, m_indexStream);
 575     elementsWritten += fwrite(&startPosition,  sizeof(startPosition),  1, m_indexStream);
 576     return ( elementsWritten == 3 );
 577 }