src/api/BamIndex.cpp

   1 // ***************************************************************************
   2 // BamIndex.cpp (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 3 September 2010 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides index functionality - both for the default (standardized) BAM
   9 // index format (.bai) as well as a BamTools-specific (nonstandard) index
  10 // format (.bti).
  11 // ***************************************************************************
  12
  13 #include <cstdio>
  14 #include <cstdlib>
  15 #include <algorithm>
  16 #include <map>
  17 #include "BamIndex.h"
  18 #include "BamReader.h"
  19 #include "BGZF.h"
  20 using namespace std;
  21 using namespace BamTools;
  22
  23 // -------------------------------
  24 // BamIndex implementation
  25
  26 BamIndex::BamIndex(BamTools::BgzfData* bgzf, BamTools::BamReader* reader, bool isBigEndian)
  27     : m_BGZF(bgzf)
  28     , m_reader(reader)
  29     , m_isBigEndian(isBigEndian)
  30 {
  31     if ( m_reader && m_reader->IsOpen() )
  32         m_references = m_reader->GetReferenceData();
  33 }
  34
  35 bool BamIndex::HasAlignments(const int& referenceID) {
  36
  37     // return false if invalid ID
  38     if ( (referenceID < 0) || (referenceID >= (int)m_references.size()) )
  39         return false;
  40
  41     // else return status of reference (has alignments?)
  42     else
  43         return m_references.at(referenceID).RefHasAlignments;
  44 }
  45
  46 // #########################################################################################
  47 // #########################################################################################
  48
  49 // -------------------------------
  50 // BamStandardIndex structs & typedefs
  51
  52 namespace BamTools {
  53
  54 // --------------------------------------------------
  55 // BamStandardIndex data structures & typedefs
  56 struct Chunk {
  57
  58     // data members
  59     uint64_t Start;
  60     uint64_t Stop;
  61
  62     // constructor
  63     Chunk(const uint64_t& start = 0,
  64           const uint64_t& stop = 0)
  65         : Start(start)
  66         , Stop(stop)
  67     { }
  68 };
  69
  70 bool ChunkLessThan(const Chunk& lhs, const Chunk& rhs) {
  71     return lhs.Start < rhs.Start;
  72 }
  73
  74 typedef vector<Chunk> ChunkVector;
  75 typedef map<uint32_t, ChunkVector> BamBinMap;
  76 typedef vector<uint64_t> LinearOffsetVector;
  77
  78 struct ReferenceIndex {
  79
  80     // data members
  81     BamBinMap Bins;
  82     LinearOffsetVector Offsets;
  83
  84     // constructor
  85     ReferenceIndex(const BamBinMap& binMap = BamBinMap(),
  86                    const LinearOffsetVector& offsets = LinearOffsetVector())
  87         : Bins(binMap)
  88         , Offsets(offsets)
  89     { }
  90 };
  91
  92 typedef vector<ReferenceIndex> BamStandardIndexData;
  93
  94 } // namespace BamTools
  95
  96 // -------------------------------
  97 // BamStandardIndex implementation
  98
  99 struct BamStandardIndex::BamStandardIndexPrivate {
 100
 101     // -------------------------
 102     // data members
 103
 104     BamStandardIndexData m_indexData;
 105     BamStandardIndex*    m_parent;
 106
 107     // -------------------------
 108     // ctor & dtor
 109
 110     BamStandardIndexPrivate(BamStandardIndex* parent) : m_parent(parent) { }
 111     ~BamStandardIndexPrivate(void) { m_indexData.clear(); }
 112
 113     // -------------------------
 114     // internal methods
 115
 116     // calculate bins that overlap region
 117     int BinsFromRegion(const BamTools::BamRegion& region, const bool isRightBoundSpecified, uint16_t bins[BamTools::MAX_BIN]);
 118     // saves BAM bin entry for index
 119     void InsertBinEntry(BamBinMap& binMap, const uint32_t& saveBin, const uint64_t& saveOffset, const uint64_t& lastOffset);
 120     // saves linear offset entry for index
 121     void InsertLinearOffset(LinearOffsetVector& offsets, const BamAlignment& bAlignment, const uint64_t& lastOffset);
 122     // simplifies index by merging 'chunks'
 123     void MergeChunks(void);
 124
 125 };
 126
 127 BamStandardIndex::BamStandardIndex(BgzfData* bgzf, BamReader* reader, bool isBigEndian)
 128     : BamIndex(bgzf, reader, isBigEndian)
 129 {
 130     d = new BamStandardIndexPrivate(this);
 131 }
 132
 133 BamStandardIndex::~BamStandardIndex(void) {
 134     delete d;
 135     d = 0;
 136 }
 137
 138 // calculate bins that overlap region
 139 int BamStandardIndex::BamStandardIndexPrivate::BinsFromRegion(const BamRegion& region, const bool isRightBoundSpecified, uint16_t bins[MAX_BIN]) {
 140
 141     // get region boundaries
 142     uint32_t begin = (unsigned int)region.LeftPosition;
 143     uint32_t end;
 144
 145     // if right bound specified AND left&right bounds are on same reference
 146     // OK to use right bound position
 147     if ( isRightBoundSpecified && ( region.LeftRefID == region.RightRefID ) )
 148         end = (unsigned int)region.RightPosition;
 149
 150     // otherwise, use end of left bound reference as cutoff
 151     else
 152         end = (unsigned int)m_parent->m_references.at(region.LeftRefID).RefLength - 1;
 153
 154     // initialize list, bin '0' always a valid bin
 155     int i = 0;
 156     bins[i++] = 0;
 157
 158     // get rest of bins that contain this region
 159     unsigned int k;
 160     for (k =    1 + (begin>>26); k <=    1 + (end>>26); ++k) { bins[i++] = k; }
 161     for (k =    9 + (begin>>23); k <=    9 + (end>>23); ++k) { bins[i++] = k; }
 162     for (k =   73 + (begin>>20); k <=   73 + (end>>20); ++k) { bins[i++] = k; }
 163     for (k =  585 + (begin>>17); k <=  585 + (end>>17); ++k) { bins[i++] = k; }
 164     for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { bins[i++] = k; }
 165
 166     // return number of bins stored
 167     return i;
 168 }
 169
 170 bool BamStandardIndex::Build(void) {
 171
 172     // be sure reader & BGZF file are valid & open for reading
 173     if ( m_reader == 0 || m_BGZF == 0 || !m_BGZF->IsOpen )
 174         return false;
 175
 176     // move file pointer to beginning of alignments
 177     m_reader->Rewind();
 178
 179     // get reference count, reserve index space
 180     int numReferences = (int)m_references.size();
 181     for ( int i = 0; i < numReferences; ++i ) {
 182         d->m_indexData.push_back(ReferenceIndex());
 183     }
 184
 185     // sets default constant for bin, ID, offset, coordinate variables
 186     const uint32_t defaultValue = 0xffffffffu;
 187
 188     // bin data
 189     uint32_t saveBin(defaultValue);
 190     uint32_t lastBin(defaultValue);
 191
 192     // reference ID data
 193     int32_t saveRefID(defaultValue);
 194     int32_t lastRefID(defaultValue);
 195
 196     // offset data
 197     uint64_t saveOffset = m_BGZF->Tell();
 198     uint64_t lastOffset = saveOffset;
 199
 200     // coordinate data
 201     int32_t lastCoordinate = defaultValue;
 202
 203     BamAlignment bAlignment;
 204     while ( m_reader->GetNextAlignmentCore(bAlignment) ) {
 205
 206         // change of chromosome, save ID, reset bin
 207         if ( lastRefID != bAlignment.RefID ) {
 208             lastRefID = bAlignment.RefID;
 209             lastBin   = defaultValue;
 210         }
 211
 212         // if lastCoordinate greater than BAM position - file not sorted properly
 213         else if ( lastCoordinate > bAlignment.Position ) {
 214             fprintf(stderr, "BAM file not properly sorted:\n");
 215             fprintf(stderr, "Alignment %s : %d > %d on reference (id = %d)", bAlignment.Name.c_str(), lastCoordinate, bAlignment.Position, bAlignment.RefID);
 216             exit(1);
 217         }
 218
 219         // if valid reference && BAM bin spans some minimum cutoff (smaller bin ids span larger regions)
 220         if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) {
 221
 222             // save linear offset entry (matched to BAM entry refID)
 223             ReferenceIndex& refIndex = d->m_indexData.at(bAlignment.RefID);
 224             LinearOffsetVector& offsets = refIndex.Offsets;
 225             d->InsertLinearOffset(offsets, bAlignment, lastOffset);
 226         }
 227
 228         // if current BamAlignment bin != lastBin, "then possibly write the binning index"
 229         if ( bAlignment.Bin != lastBin ) {
 230
 231             // if not first time through
 232             if ( saveBin != defaultValue ) {
 233
 234                 // save Bam bin entry
 235                 ReferenceIndex& refIndex = d->m_indexData.at(saveRefID);
 236                 BamBinMap& binMap = refIndex.Bins;
 237                 d->InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);
 238             }
 239
 240             // update saveOffset
 241             saveOffset = lastOffset;
 242
 243             // update bin values
 244             saveBin = bAlignment.Bin;
 245             lastBin = bAlignment.Bin;
 246
 247             // update saveRefID
 248             saveRefID = bAlignment.RefID;
 249
 250             // if invalid RefID, break out
 251             if ( saveRefID < 0 ) break;
 252         }
 253
 254         // make sure that current file pointer is beyond lastOffset
 255         if ( m_BGZF->Tell() <= (int64_t)lastOffset ) {
 256             fprintf(stderr, "Error in BGZF offsets.\n");
 257             exit(1);
 258         }
 259
 260         // update lastOffset
 261         lastOffset = m_BGZF->Tell();
 262
 263         // update lastCoordinate
 264         lastCoordinate = bAlignment.Position;
 265     }
 266
 267     // save any leftover BAM data (as long as refID is valid)
 268     if ( saveRefID >= 0 ) {
 269         // save Bam bin entry
 270         ReferenceIndex& refIndex = d->m_indexData.at(saveRefID);
 271         BamBinMap& binMap = refIndex.Bins;
 272         d->InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);
 273     }
 274
 275     // simplify index by merging chunks
 276     d->MergeChunks();
 277
 278     // iterate through references in index
 279     // store whether reference has data &
 280     // sort offsets in linear offset vector
 281     BamStandardIndexData::iterator indexIter = d->m_indexData.begin();
 282     BamStandardIndexData::iterator indexEnd  = d->m_indexData.end();
 283     for ( int i = 0; indexIter != indexEnd; ++indexIter, ++i ) {
 284
 285         // get reference index data
 286         ReferenceIndex& refIndex = (*indexIter);
 287         BamBinMap& binMap = refIndex.Bins;
 288         LinearOffsetVector& offsets = refIndex.Offsets;
 289
 290         // store whether reference has alignments or no
 291         m_references[i].RefHasAlignments = ( binMap.size() > 0 );
 292
 293         // sort linear offsets
 294         sort(offsets.begin(), offsets.end());
 295     }
 296
 297     // rewind file pointer to beginning of alignments, return success/fail
 298     return m_reader->Rewind();
 299 }
 300
 301 bool BamStandardIndex::GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets) {
 302
 303     // calculate which bins overlap this region
 304     uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);
 305     int numBins = d->BinsFromRegion(region, isRightBoundSpecified, bins);
 306
 307     // get bins for this reference
 308     const ReferenceIndex& refIndex = d->m_indexData.at(region.LeftRefID);
 309     const BamBinMap& binMap        = refIndex.Bins;
 310
 311     // get minimum offset to consider
 312     const LinearOffsetVector& linearOffsets = refIndex.Offsets;
 313     uint64_t minOffset = ( (unsigned int)(region.LeftPosition>>BAM_LIDX_SHIFT) >= linearOffsets.size() ) ? 0 : linearOffsets.at(region.LeftPosition>>BAM_LIDX_SHIFT);
 314
 315     // store all alignment 'chunk' starts (file offsets) for bins in this region
 316     for ( int i = 0; i < numBins; ++i ) {
 317
 318         const uint16_t binKey = bins[i];
 319         map<uint32_t, ChunkVector>::const_iterator binIter = binMap.find(binKey);
 320         if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) {
 321
 322             // iterate over chunks
 323             const ChunkVector& chunks = (*binIter).second;
 324             std::vector<Chunk>::const_iterator chunksIter = chunks.begin();
 325             std::vector<Chunk>::const_iterator chunksEnd  = chunks.end();
 326             for ( ; chunksIter != chunksEnd; ++chunksIter) {
 327
 328                 // if valid chunk found, store its file offset
 329                 const Chunk& chunk = (*chunksIter);
 330                 if ( chunk.Stop > minOffset )
 331                     offsets.push_back( chunk.Start );
 332             }
 333         }
 334     }
 335
 336     // clean up memory
 337     free(bins);
 338
 339     // sort the offsets before returning
 340     sort(offsets.begin(), offsets.end());
 341
 342     // return whether any offsets were found
 343     return ( offsets.size() != 0 );
 344 }
 345
 346 // saves BAM bin entry for index
 347 void BamStandardIndex::BamStandardIndexPrivate::InsertBinEntry(BamBinMap&      binMap,
 348                                      const uint32_t& saveBin,
 349                                      const uint64_t& saveOffset,
 350                                      const uint64_t& lastOffset)
 351 {
 352     // look up saveBin
 353     BamBinMap::iterator binIter = binMap.find(saveBin);
 354
 355     // create new chunk
 356     Chunk newChunk(saveOffset, lastOffset);
 357
 358     // if entry doesn't exist
 359     if ( binIter == binMap.end() ) {
 360         ChunkVector newChunks;
 361         newChunks.push_back(newChunk);
 362         binMap.insert( pair<uint32_t, ChunkVector>(saveBin, newChunks));
 363     }
 364
 365     // otherwise
 366     else {
 367         ChunkVector& binChunks = (*binIter).second;
 368         binChunks.push_back( newChunk );
 369     }
 370 }
 371
 372 // saves linear offset entry for index
 373 void BamStandardIndex::BamStandardIndexPrivate::InsertLinearOffset(LinearOffsetVector& offsets,
 374                                         const BamAlignment& bAlignment,
 375                                         const uint64_t&     lastOffset)
 376 {
 377     // get converted offsets
 378     int beginOffset = bAlignment.Position >> BAM_LIDX_SHIFT;
 379     int endOffset   = (bAlignment.GetEndPosition() - 1) >> BAM_LIDX_SHIFT;
 380
 381     // resize vector if necessary
 382     int oldSize = offsets.size();
 383     int newSize = endOffset + 1;
 384     if ( oldSize < newSize )
 385         offsets.resize(newSize, 0);
 386
 387     // store offset
 388     for( int i = beginOffset + 1; i <= endOffset; ++i ) {
 389         if ( offsets[i] == 0 )
 390             offsets[i] = lastOffset;
 391     }
 392 }
 393
 394 bool BamStandardIndex::Load(const string& filename)  {
 395
 396     // open index file, abort on error
 397     FILE* indexStream = fopen(filename.c_str(), "rb");
 398     if( !indexStream ) {
 399         fprintf(stderr, "ERROR: Unable to open the BAM index file %s for reading.\n", filename.c_str());
 400         return false;
 401     }
 402
 403     // set placeholder to receive input byte count (suppresses compiler warnings)
 404     size_t elementsRead = 0;
 405
 406     // see if index is valid BAM index
 407     char magic[4];
 408     elementsRead = fread(magic, 1, 4, indexStream);
 409     if ( strncmp(magic, "BAI\1", 4) ) {
 410         fprintf(stderr, "Problem with index file - invalid format.\n");
 411         fclose(indexStream);
 412         return false;
 413     }
 414
 415     // get number of reference sequences
 416     uint32_t numRefSeqs;
 417     elementsRead = fread(&numRefSeqs, 4, 1, indexStream);
 418     if ( m_isBigEndian ) SwapEndian_32(numRefSeqs);
 419
 420     // intialize space for BamStandardIndexData data structure
 421     d->m_indexData.reserve(numRefSeqs);
 422
 423     // iterate over reference sequences
 424     for ( unsigned int i = 0; i < numRefSeqs; ++i ) {
 425
 426         // get number of bins for this reference sequence
 427         int32_t numBins;
 428         elementsRead = fread(&numBins, 4, 1, indexStream);
 429         if ( m_isBigEndian ) SwapEndian_32(numBins);
 430
 431         if ( numBins > 0 ) {
 432             RefData& refEntry = m_references[i];
 433             refEntry.RefHasAlignments = true;
 434         }
 435
 436         // intialize BinVector
 437         BamBinMap binMap;
 438
 439         // iterate over bins for that reference sequence
 440         for ( int j = 0; j < numBins; ++j ) {
 441
 442             // get binID
 443             uint32_t binID;
 444             elementsRead = fread(&binID, 4, 1, indexStream);
 445
 446             // get number of regionChunks in this bin
 447             uint32_t numChunks;
 448             elementsRead = fread(&numChunks, 4, 1, indexStream);
 449
 450             if ( m_isBigEndian ) {
 451               SwapEndian_32(binID);
 452               SwapEndian_32(numChunks);
 453             }
 454
 455             // intialize ChunkVector
 456             ChunkVector regionChunks;
 457             regionChunks.reserve(numChunks);
 458
 459             // iterate over regionChunks in this bin
 460             for ( unsigned int k = 0; k < numChunks; ++k ) {
 461
 462                 // get chunk boundaries (left, right)
 463                 uint64_t left;
 464                 uint64_t right;
 465                 elementsRead = fread(&left, 8, 1, indexStream);
 466                 elementsRead = fread(&right, 8, 1, indexStream);
 467
 468                 if ( m_isBigEndian ) {
 469                     SwapEndian_64(left);
 470                     SwapEndian_64(right);
 471                 }
 472
 473                 // save ChunkPair
 474                 regionChunks.push_back( Chunk(left, right) );
 475             }
 476
 477             // sort chunks for this bin
 478             sort( regionChunks.begin(), regionChunks.end(), ChunkLessThan );
 479
 480             // save binID, chunkVector for this bin
 481             binMap.insert( pair<uint32_t, ChunkVector>(binID, regionChunks) );
 482         }
 483
 484         // -----------------------------------------------------
 485         // load linear index for this reference sequence
 486
 487         // get number of linear offsets
 488         int32_t numLinearOffsets;
 489         elementsRead = fread(&numLinearOffsets, 4, 1, indexStream);
 490         if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets);
 491
 492         // intialize LinearOffsetVector
 493         LinearOffsetVector offsets;
 494         offsets.reserve(numLinearOffsets);
 495
 496         // iterate over linear offsets for this reference sequeence
 497         uint64_t linearOffset;
 498         for ( int j = 0; j < numLinearOffsets; ++j ) {
 499             // read a linear offset & store
 500             elementsRead = fread(&linearOffset, 8, 1, indexStream);
 501             if ( m_isBigEndian ) SwapEndian_64(linearOffset);
 502             offsets.push_back(linearOffset);
 503         }
 504
 505         // sort linear offsets
 506         sort( offsets.begin(), offsets.end() );
 507
 508         // store index data for that reference sequence
 509         d->m_indexData.push_back( ReferenceIndex(binMap, offsets) );
 510     }
 511
 512     // close index file (.bai) and return
 513     fclose(indexStream);
 514     return true;
 515 }
 516
 517 // merges 'alignment chunks' in BAM bin (used for index building)
 518 void BamStandardIndex::BamStandardIndexPrivate::MergeChunks(void) {
 519
 520     // iterate over reference enties
 521     BamStandardIndexData::iterator indexIter = m_indexData.begin();
 522     BamStandardIndexData::iterator indexEnd  = m_indexData.end();
 523     for ( ; indexIter != indexEnd; ++indexIter ) {
 524
 525         // get BAM bin map for this reference
 526         ReferenceIndex& refIndex = (*indexIter);
 527         BamBinMap& bamBinMap = refIndex.Bins;
 528
 529         // iterate over BAM bins
 530         BamBinMap::iterator binIter = bamBinMap.begin();
 531         BamBinMap::iterator binEnd  = bamBinMap.end();
 532         for ( ; binIter != binEnd; ++binIter ) {
 533
 534             // get chunk vector for this bin
 535             ChunkVector& binChunks = (*binIter).second;
 536             if ( binChunks.size() == 0 ) continue;
 537
 538             ChunkVector mergedChunks;
 539             mergedChunks.push_back( binChunks[0] );
 540
 541             // iterate over chunks
 542             int i = 0;
 543             ChunkVector::iterator chunkIter = binChunks.begin();
 544             ChunkVector::iterator chunkEnd  = binChunks.end();
 545             for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {
 546
 547                 // get 'currentChunk' based on numeric index
 548                 Chunk& currentChunk = mergedChunks[i];
 549
 550                 // get iteratorChunk based on vector iterator
 551                 Chunk& iteratorChunk = (*chunkIter);
 552
 553                 // if currentChunk.Stop(shifted) == iterator Chunk.Start(shifted)
 554                 if ( currentChunk.Stop>>16 == iteratorChunk.Start>>16 ) {
 555
 556                     // set currentChunk.Stop to iteratorChunk.Stop
 557                     currentChunk.Stop = iteratorChunk.Stop;
 558                 }
 559
 560                 // otherwise
 561                 else {
 562                     // set currentChunk + 1 to iteratorChunk
 563                     mergedChunks.push_back(iteratorChunk);
 564                     ++i;
 565                 }
 566             }
 567
 568             // saved merged chunk vector
 569             (*binIter).second = mergedChunks;
 570         }
 571     }
 572 }
 573
 574 // writes in-memory index data out to file
 575 // N.B. - (this is the original BAM filename, method will modify it to use applicable extension)
 576 bool BamStandardIndex::Write(const std::string& bamFilename) {
 577
 578     string indexFilename = bamFilename + ".bai";
 579     FILE* indexStream = fopen(indexFilename.c_str(), "wb");
 580     if ( indexStream == 0 ) {
 581         fprintf(stderr, "ERROR: Could not open file to save index.\n");
 582         return false;
 583     }
 584
 585     // write BAM index header
 586     fwrite("BAI\1", 1, 4, indexStream);
 587
 588     // write number of reference sequences
 589     int32_t numReferenceSeqs = d->m_indexData.size();
 590     if ( m_isBigEndian ) SwapEndian_32(numReferenceSeqs);
 591     fwrite(&numReferenceSeqs, 4, 1, indexStream);
 592
 593     // iterate over reference sequences
 594     BamStandardIndexData::const_iterator indexIter = d->m_indexData.begin();
 595     BamStandardIndexData::const_iterator indexEnd  = d->m_indexData.end();
 596     for ( ; indexIter != indexEnd; ++ indexIter ) {
 597
 598         // get reference index data
 599         const ReferenceIndex& refIndex = (*indexIter);
 600         const BamBinMap& binMap = refIndex.Bins;
 601         const LinearOffsetVector& offsets = refIndex.Offsets;
 602
 603         // write number of bins
 604         int32_t binCount = binMap.size();
 605         if ( m_isBigEndian ) SwapEndian_32(binCount);
 606         fwrite(&binCount, 4, 1, indexStream);
 607
 608         // iterate over bins
 609         BamBinMap::const_iterator binIter = binMap.begin();
 610         BamBinMap::const_iterator binEnd  = binMap.end();
 611         for ( ; binIter != binEnd; ++binIter ) {
 612
 613             // get bin data (key and chunk vector)
 614             uint32_t binKey = (*binIter).first;
 615             const ChunkVector& binChunks = (*binIter).second;
 616
 617             // save BAM bin key
 618             if ( m_isBigEndian ) SwapEndian_32(binKey);
 619             fwrite(&binKey, 4, 1, indexStream);
 620
 621             // save chunk count
 622             int32_t chunkCount = binChunks.size();
 623             if ( m_isBigEndian ) SwapEndian_32(chunkCount);
 624             fwrite(&chunkCount, 4, 1, indexStream);
 625
 626             // iterate over chunks
 627             ChunkVector::const_iterator chunkIter = binChunks.begin();
 628             ChunkVector::const_iterator chunkEnd  = binChunks.end();
 629             for ( ; chunkIter != chunkEnd; ++chunkIter ) {
 630
 631                 // get current chunk data
 632                 const Chunk& chunk    = (*chunkIter);
 633                 uint64_t start = chunk.Start;
 634                 uint64_t stop  = chunk.Stop;
 635
 636                 if ( m_isBigEndian ) {
 637                     SwapEndian_64(start);
 638                     SwapEndian_64(stop);
 639                 }
 640
 641                 // save chunk offsets
 642                 fwrite(&start, 8, 1, indexStream);
 643                 fwrite(&stop,  8, 1, indexStream);
 644             }
 645         }
 646
 647         // write linear offsets size
 648         int32_t offsetSize = offsets.size();
 649         if ( m_isBigEndian ) SwapEndian_32(offsetSize);
 650         fwrite(&offsetSize, 4, 1, indexStream);
 651
 652         // iterate over linear offsets
 653         LinearOffsetVector::const_iterator offsetIter = offsets.begin();
 654         LinearOffsetVector::const_iterator offsetEnd  = offsets.end();
 655         for ( ; offsetIter != offsetEnd; ++offsetIter ) {
 656
 657             // write linear offset value
 658             uint64_t linearOffset = (*offsetIter);
 659             if ( m_isBigEndian ) SwapEndian_64(linearOffset);
 660             fwrite(&linearOffset, 8, 1, indexStream);
 661         }
 662     }
 663
 664     // flush buffer, close file, and return success
 665     fflush(indexStream);
 666     fclose(indexStream);
 667     return true;
 668 }
 669
 670 // #########################################################################################
 671 // #########################################################################################
 672
 673 // -------------------------------------
 674 // BamToolsIndex implementation
 675
 676 namespace BamTools {
 677
 678 struct BamToolsIndexEntry {
 679
 680     // data members
 681     int64_t Offset;
 682     int RefID;
 683     int Position;
 684
 685     // ctor
 686     BamToolsIndexEntry(const uint64_t& offset = 0,
 687                        const int& id = -1,
 688                        const int& position = -1)
 689         : Offset(offset)
 690         , RefID(id)
 691         , Position(position)
 692     { }
 693 };
 694
 695 typedef vector<BamToolsIndexEntry> BamToolsIndexData;
 696
 697 } // namespace BamTools
 698
 699 struct BamToolsIndex::BamToolsIndexPrivate {
 700
 701     // -------------------------
 702     // data members
 703     BamToolsIndexData m_indexData;
 704     BamToolsIndex*    m_parent;
 705     int32_t           m_blockSize;
 706
 707     // -------------------------
 708     // ctor & dtor
 709
 710     BamToolsIndexPrivate(BamToolsIndex* parent)
 711         : m_parent(parent)
 712         , m_blockSize(1000)
 713     { }
 714
 715     ~BamToolsIndexPrivate(void) { }
 716
 717     // -------------------------
 718     // internal methods
 719 };
 720
 721 BamToolsIndex::BamToolsIndex(BgzfData* bgzf, BamReader* reader, bool isBigEndian)
 722     : BamIndex(bgzf, reader, isBigEndian)
 723 {
 724     d = new BamToolsIndexPrivate(this);
 725 }
 726
 727 BamToolsIndex::~BamToolsIndex(void) {
 728     delete d;
 729     d = 0;
 730 }
 731
 732 bool BamToolsIndex::Build(void) {
 733
 734     // be sure reader & BGZF file are valid & open for reading
 735     if ( m_reader == 0 || m_BGZF == 0 || !m_BGZF->IsOpen )
 736         return false;
 737
 738     // move file pointer to beginning of alignments
 739     m_reader->Rewind();
 740
 741     // plow through alignments, store block offsets
 742     int32_t currentBlockCount  = 0;
 743     int64_t blockStartOffset   = m_BGZF->Tell();
 744     int     blockStartId       = -1;
 745     int     blockStartPosition = -1;
 746     BamAlignment al;
 747     while ( m_reader->GetNextAlignmentCore(al) ) {
 748
 749         // set reference flag
 750         m_references[al.RefID].RefHasAlignments = true;
 751
 752         // if beginning of block, save first alignment's refID & position
 753         if ( currentBlockCount == 0 ) {
 754             blockStartId = al.RefID;
 755             blockStartPosition = al.Position;
 756         }
 757
 758         // increment block counter
 759         ++currentBlockCount;
 760
 761         // if block is full, get offset for next block, reset currentBlockCount
 762         if ( currentBlockCount == d->m_blockSize ) {
 763             d->m_indexData.push_back( BamToolsIndexEntry(blockStartOffset, blockStartId, blockStartPosition) );
 764             blockStartOffset = m_BGZF->Tell();
 765             currentBlockCount = 0;
 766         }
 767     }
 768
 769     return m_reader->Rewind();
 770 }
 771
 772 // N.B. - ignores isRightBoundSpecified
 773 bool BamToolsIndex::GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets) {
 774
 775     // return false if no index data present
 776     if ( d->m_indexData.empty() ) return false;
 777
 778     // clear any prior data
 779     offsets.clear();
 780
 781     // calculate nearest index to jump to
 782     int64_t previousOffset = -1;
 783     BamToolsIndexData::const_iterator indexIter = d->m_indexData.begin();
 784     BamToolsIndexData::const_iterator indexEnd  = d->m_indexData.end();
 785     for ( ; indexIter != indexEnd; ++indexIter ) {
 786
 787         const BamToolsIndexEntry& entry = (*indexIter);
 788
 789         // check if we are 'past' beginning of desired region
 790         // if so, we will break out & use previously stored offset
 791         if ( entry.RefID > region.LeftRefID ) break;
 792         if ( (entry.RefID == region.LeftRefID) && (entry.Position > region.LeftPosition) ) break;
 793
 794         // not past desired region, so store current entry offset in previousOffset
 795         previousOffset = entry.Offset;
 796     }
 797
 798     // no index was found
 799     if ( previousOffset == -1 ) return false;
 800
 801     // store offset & return success
 802     offsets.push_back(previousOffset);
 803     return true;
 804 }
 805
 806 bool BamToolsIndex::Load(const string& filename) {
 807
 808     // open index file, abort on error
 809     FILE* indexStream = fopen(filename.c_str(), "rb");
 810     if( !indexStream ) {
 811         fprintf(stderr, "ERROR: Unable to open the BAM index file %s for reading.\n", filename.c_str());
 812         return false;
 813     }
 814
 815     // set placeholder to receive input byte count (suppresses compiler warnings)
 816     size_t elementsRead = 0;
 817
 818     // see if index is valid BAM index
 819     char magic[4];
 820     elementsRead = fread(magic, 1, 4, indexStream);
 821     if ( strncmp(magic, "BTI\1", 4) ) {
 822         fprintf(stderr, "Problem with index file - invalid format.\n");
 823         fclose(indexStream);
 824         return false;
 825     }
 826
 827     // read in block size
 828     elementsRead = fread(&d->m_blockSize, sizeof(d->m_blockSize), 1, indexStream);
 829     if ( m_isBigEndian ) SwapEndian_32(d->m_blockSize);
 830
 831     // read in number of offsets
 832     uint32_t numOffsets;
 833     elementsRead = fread(&numOffsets, sizeof(numOffsets), 1, indexStream);
 834     if ( m_isBigEndian ) SwapEndian_32(numOffsets);
 835
 836     // reserve space for index data
 837     d->m_indexData.reserve(numOffsets);
 838
 839     // iterate over index entries
 840     for ( unsigned int i = 0; i < numOffsets; ++i ) {
 841
 842         uint64_t offset;
 843         int id;
 844         int position;
 845
 846         // read in data
 847         elementsRead = fread(&offset, sizeof(offset), 1, indexStream);
 848         elementsRead = fread(&id, sizeof(id), 1, indexStream);
 849         elementsRead = fread(&position, sizeof(position), 1, indexStream);
 850
 851         // swap endian-ness if necessary
 852         if ( m_isBigEndian ) {
 853             SwapEndian_64(offset);
 854             SwapEndian_32(id);
 855             SwapEndian_32(position);
 856         }
 857
 858         // save reference index entry
 859         d->m_indexData.push_back( BamToolsIndexEntry(offset, id, position) );
 860
 861         // set reference flag
 862         m_references[id].RefHasAlignments = true;       // what about sparse references? wont be able to set flag?
 863     }
 864
 865     // close index file and return
 866     fclose(indexStream);
 867     return true;
 868 }
 869
 870 // writes in-memory index data out to file
 871 // N.B. - (this is the original BAM filename, method will modify it to use applicable extension)
 872 bool BamToolsIndex::Write(const std::string& bamFilename) {
 873
 874     string indexFilename = bamFilename + ".bti";
 875     FILE* indexStream = fopen(indexFilename.c_str(), "wb");
 876     if ( indexStream == 0 ) {
 877         fprintf(stderr, "ERROR: Could not open file to save index.\n");
 878         return false;
 879     }
 880
 881     // write BAM index header
 882     fwrite("BTI\1", 1, 4, indexStream);
 883
 884     // write block size
 885     int32_t blockSize = d->m_blockSize;
 886     if ( m_isBigEndian ) SwapEndian_32(blockSize);
 887     fwrite(&blockSize, sizeof(blockSize), 1, indexStream);
 888
 889     // write number of offset entries
 890     uint32_t numOffsets = d->m_indexData.size();
 891     if ( m_isBigEndian ) SwapEndian_32(numOffsets);
 892     fwrite(&numOffsets, sizeof(numOffsets), 1, indexStream);
 893
 894     // iterate over offset entries
 895     BamToolsIndexData::const_iterator indexIter = d->m_indexData.begin();
 896     BamToolsIndexData::const_iterator indexEnd  = d->m_indexData.end();
 897     for ( ; indexIter != indexEnd; ++ indexIter ) {
 898
 899         // get reference index data
 900         const BamToolsIndexEntry& entry = (*indexIter);
 901
 902         // copy entry data
 903         uint64_t offset = entry.Offset;
 904         int id = entry.RefID;
 905         int position = entry.Position;
 906
 907         // swap endian-ness if necessary
 908         if ( m_isBigEndian ) {
 909             SwapEndian_64(offset);
 910             SwapEndian_32(id);
 911             SwapEndian_32(position);
 912         }
 913
 914         // write the reference index entry
 915         fwrite(&offset,   sizeof(offset), 1, indexStream);
 916         fwrite(&id,       sizeof(id), 1, indexStream);
 917         fwrite(&position, sizeof(position), 1, indexStream);
 918     }
 919
 920     // flush file buffer, close file, and return success
 921     fflush(indexStream);
 922     fclose(indexStream);
 923     return true;
 924 }