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