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