src/api/internal/BamStandardIndex_p.cpp

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