src/api/internal/BamStandardIndex_p.cpp

   1 // ***************************************************************************
   2 // BamStandardIndex.cpp (c) 2010 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 24 June 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides index operations for the standardized BAM index format (".bai")
   8 // ***************************************************************************
   9
  10 #include <api/BamAlignment.h>
  11 #include <api/internal/BamReader_p.h>
  12 #include <api/internal/BamStandardIndex_p.h>
  13 using namespace BamTools;
  14 using namespace BamTools::Internal;
  15
  16 #include <cstdio>
  17 #include <cstdlib>
  18 #include <cstring>
  19 #include <algorithm>
  20 #include <iostream>
  21 using namespace std;
  22
  23 // static BamStandardIndex constants
  24 const int BamStandardIndex::MAX_BIN               = 37450;  // =(8^6-1)/7+1
  25 const int BamStandardIndex::BAM_LIDX_SHIFT        = 14;
  26 const string BamStandardIndex::BAI_EXTENSION      = ".bai";
  27 const char* const BamStandardIndex::BAI_MAGIC     = "BAI\1";
  28 const int BamStandardIndex::SIZEOF_ALIGNMENTCHUNK = sizeof(uint64_t)*2;
  29 const int BamStandardIndex::SIZEOF_BINCORE        = sizeof(uint32_t) + sizeof(int32_t);
  30 const int BamStandardIndex::SIZEOF_LINEAROFFSET   = sizeof(uint64_t);
  31
  32 // ctor
  33 BamStandardIndex::BamStandardIndex(Internal::BamReaderPrivate* reader)
  34     : BamIndex(reader)
  35     , m_indexStream(0)
  36     , m_cacheMode(BamIndex::LimitedIndexCaching)
  37     , m_buffer(0)
  38     , m_bufferLength(0)
  39 {
  40      m_isBigEndian = BamTools::SystemIsBigEndian();
  41 }
  42
  43 // dtor
  44 BamStandardIndex::~BamStandardIndex(void) {
  45     CloseFile();
  46 }
  47
  48 bool BamStandardIndex::AdjustRegion(const BamRegion& region, uint32_t& begin, uint32_t& end) {
  49
  50     // retrieve references from reader
  51     const RefVector& references = m_reader->GetReferenceData();
  52
  53     // make sure left-bound position is valid
  54     if ( region.LeftPosition > references.at(region.LeftRefID).RefLength )
  55         return false;
  56
  57     // set region 'begin'
  58     begin = (unsigned int)region.LeftPosition;
  59
  60     // if right bound specified AND left&right bounds are on same reference
  61     // OK to use right bound position as region 'end'
  62     if ( region.isRightBoundSpecified() && ( region.LeftRefID == region.RightRefID ) )
  63         end = (unsigned int)region.RightPosition;
  64
  65     // otherwise, set region 'end' to last reference base
  66     else end = (unsigned int)references.at(region.LeftRefID).RefLength - 1;
  67
  68     // return success
  69     return true;
  70 }
  71
  72 void BamStandardIndex::CalculateCandidateBins(const uint32_t& begin,
  73                                               const uint32_t& end,
  74                                               set<uint16_t>& candidateBins)
  75 {
  76     // initialize list, bin '0' is always a valid bin
  77     candidateBins.insert(0);
  78
  79     // get rest of bins that contain this region
  80     unsigned int k;
  81     for (k =    1 + (begin>>26); k <=    1 + (end>>26); ++k) { candidateBins.insert(k); }
  82     for (k =    9 + (begin>>23); k <=    9 + (end>>23); ++k) { candidateBins.insert(k); }
  83     for (k =   73 + (begin>>20); k <=   73 + (end>>20); ++k) { candidateBins.insert(k); }
  84     for (k =  585 + (begin>>17); k <=  585 + (end>>17); ++k) { candidateBins.insert(k); }
  85     for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { candidateBins.insert(k); }
  86 }
  87
  88 bool BamStandardIndex::CalculateCandidateOffsets(const BaiReferenceSummary& refSummary,
  89                                                  const uint64_t& minOffset,
  90                                                  set<uint16_t>& candidateBins,
  91                                                  vector<int64_t>& offsets)
  92 {
  93     // attempt seek to first bin
  94     if ( !Seek(refSummary.FirstBinFilePosition, SEEK_SET) )
  95         return false;
  96
  97     // iterate over reference bins
  98     uint32_t binId;
  99     int32_t numAlignmentChunks;
 100     set<uint16_t>::iterator candidateBinIter;
 101     for ( int i = 0; i < refSummary.NumBins; ++i ) {
 102
 103         // read bin contents (if successful, alignment chunks are now in m_buffer)
 104         if ( !ReadBinIntoBuffer(binId, numAlignmentChunks) )
 105             return false;
 106
 107         // see if bin is a 'candidate bin'
 108         candidateBinIter = candidateBins.find(binId);
 109
 110         // if not, move on to next bin
 111         if ( candidateBinIter == candidateBins.end() )
 112             continue;
 113
 114         // otherwise, check bin's contents against for overlap
 115         else {
 116
 117             unsigned int offset = 0;
 118             uint64_t chunkStart;
 119             uint64_t chunkStop;
 120
 121             // iterate over alignment chunks
 122             for (int j = 0; j < numAlignmentChunks; ++j ) {
 123
 124                 // read chunk start & stop from buffer
 125                 memcpy((char*)&chunkStart, m_buffer+offset, sizeof(uint64_t));
 126                 offset += sizeof(uint64_t);
 127                 memcpy((char*)&chunkStop, m_buffer+offset, sizeof(uint64_t));
 128                 offset += sizeof(uint64_t);
 129
 130                 // swap endian-ness if necessary
 131                 if ( m_isBigEndian ) {
 132                     SwapEndian_64(chunkStart);
 133                     SwapEndian_64(chunkStop);
 134                 }
 135
 136                 // store alignment chunk's start offset
 137                 // if its stop offset is larger than our 'minOffset'
 138                 if ( chunkStop >= minOffset )
 139                     offsets.push_back(chunkStart);
 140             }
 141
 142             // 'pop' bin ID from candidate bins set
 143             candidateBins.erase(candidateBinIter);
 144
 145             // quit if no more candidates
 146             if ( candidateBins.empty() )
 147                 break;
 148         }
 149     }
 150
 151     // return success
 152     return true;
 153 }
 154
 155 uint64_t BamStandardIndex::CalculateMinOffset(const BaiReferenceSummary& refSummary,
 156                                               const uint32_t& begin)
 157 {
 158     // if no linear offsets exist, return 0
 159     if ( refSummary.NumLinearOffsets == 0 )
 160         return 0;
 161
 162     // if 'begin' starts beyond last linear offset, use the last linear offset as minimum
 163     // else use the offset corresponding to the requested start position
 164     const int shiftedBegin = begin>>BamStandardIndex::BAM_LIDX_SHIFT;
 165     if ( shiftedBegin >= refSummary.NumLinearOffsets )
 166         return LookupLinearOffset( refSummary, refSummary.NumLinearOffsets-1 );
 167     else
 168         return LookupLinearOffset( refSummary, shiftedBegin );
 169 }
 170
 171 void BamStandardIndex::CheckBufferSize(char*& buffer,
 172                                        unsigned int& bufferLength,
 173                                        const unsigned int& requestedBytes)
 174 {
 175     try {
 176         if ( requestedBytes > bufferLength ) {
 177             bufferLength = requestedBytes + 10;
 178             delete[] buffer;
 179             buffer = new char[bufferLength];
 180         }
 181     } catch ( std::bad_alloc ) {
 182         cerr << "BamStandardIndex ERROR: out of memory when allocating "
 183              << requestedBytes << " byes" << endl;
 184         exit(1);
 185     }
 186 }
 187
 188 void BamStandardIndex::CheckBufferSize(unsigned char*& buffer,
 189                                        unsigned int& bufferLength,
 190                                        const unsigned int& requestedBytes)
 191 {
 192     try {
 193         if ( requestedBytes > bufferLength ) {
 194             bufferLength = requestedBytes + 10;
 195             delete[] buffer;
 196             buffer = new unsigned char[bufferLength];
 197         }
 198     } catch ( std::bad_alloc ) {
 199         cerr << "BamStandardIndex ERROR: out of memory when allocating "
 200              << requestedBytes << " byes" << endl;
 201         exit(1);
 202     }
 203 }
 204
 205 bool BamStandardIndex::CheckMagicNumber(void) {
 206
 207     // check 'magic number' to see if file is BAI index
 208     char magic[4];
 209     size_t elementsRead = fread(magic, sizeof(char), 4, m_indexStream);
 210     if ( elementsRead != 4 ) {
 211         cerr << "BamStandardIndex ERROR: could not read format 'magic number'" << endl;
 212         return false;
 213     }
 214
 215     // compare to expected value
 216     if ( strncmp(magic, BamStandardIndex::BAI_MAGIC, 4) != 0 ) {
 217         cerr << "BamStandardIndex ERROR: invalid format" << endl;
 218         return false;
 219     }
 220
 221     // otherwise OK
 222     return true;
 223 }
 224
 225 void BamStandardIndex::ClearReferenceEntry(BaiReferenceEntry& refEntry) {
 226     refEntry.ID = -1;
 227     refEntry.Bins.clear();
 228     refEntry.LinearOffsets.clear();
 229 }
 230
 231 void BamStandardIndex::CloseFile(void) {
 232
 233     // close file stream
 234     if ( IsFileOpen() )
 235         fclose(m_indexStream);
 236
 237     // clear index file summary data
 238     m_indexFileSummary.clear();
 239
 240     // clean up I/O buffer
 241     delete[] m_buffer;
 242     m_buffer = 0;
 243     m_bufferLength = 0;
 244 }
 245
 246 // builds index from associated BAM file & writes out to index file
 247 bool BamStandardIndex::Create(void) {
 248
 249     // return false if BamReader is invalid or not open
 250     if ( m_reader == 0 || !m_reader->IsOpen() ) {
 251         cerr << "BamStandardIndex ERROR: BamReader is not open"
 252              << ", aborting index creation" << endl;
 253         return false;
 254     }
 255
 256     // rewind BamReader
 257     if ( !m_reader->Rewind() ) {
 258         cerr << "BamStandardIndex ERROR: could not rewind BamReader to create index"
 259              << ", aborting index creation" << endl;
 260         return false;
 261     }
 262
 263     // open new index file (read & write)
 264     string indexFilename = m_reader->Filename() + Extension();
 265     if ( !OpenFile(indexFilename, "w+b") ) {
 266         cerr << "BamStandardIndex ERROR: could not open output index file: " << indexFilename
 267              << ", aborting index creation" << endl;
 268         return false;
 269     }
 270
 271     // initialize BaiFileSummary with number of references
 272     const int& numReferences = m_reader->GetReferenceCount();
 273     ReserveForSummary(numReferences);
 274
 275     // initialize output file
 276     bool createdOk = true;
 277     createdOk &= WriteHeader();
 278
 279     // set up bin, ID, offset, & coordinate markers
 280     const uint32_t defaultValue = 0xffffffffu;
 281     uint32_t currentBin    = defaultValue;
 282     uint32_t lastBin       = defaultValue;
 283     int32_t  currentRefID  = defaultValue;
 284     int32_t  lastRefID     = defaultValue;
 285     uint64_t currentOffset = (uint64_t)m_reader->Tell();
 286     uint64_t lastOffset    = currentOffset;
 287     int32_t  lastPosition  = defaultValue;
 288
 289     // iterate through alignments in BAM file
 290     BamAlignment al;
 291     BaiReferenceEntry refEntry;
 292     while ( m_reader->LoadNextAlignment(al) ) {
 293
 294         // changed to new reference
 295         if ( lastRefID != al.RefID ) {
 296
 297             // if not first reference, save previous reference data
 298             if ( lastRefID != (int32_t)defaultValue ) {
 299
 300                 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
 301                 createdOk &= WriteReferenceEntry(refEntry);
 302                 ClearReferenceEntry(refEntry);
 303
 304                 // write any empty references between (but *NOT* including) lastRefID & al.RefID
 305                 for ( int i = lastRefID+1; i < al.RefID; ++i ) {
 306                     BaiReferenceEntry emptyEntry(i);
 307                     createdOk &= WriteReferenceEntry(emptyEntry);
 308                 }
 309
 310                 // update bin markers
 311                 currentOffset = lastOffset;
 312                 currentBin    = al.Bin;
 313                 lastBin       = al.Bin;
 314                 currentRefID  = al.RefID;
 315             }
 316
 317             // first pass
 318             // write any empty references up to (but *NOT* including) al.RefID
 319             else {
 320                 for ( int i = 0; i < al.RefID; ++i ) {
 321                     BaiReferenceEntry emptyEntry(i);
 322                     createdOk &= WriteReferenceEntry(emptyEntry);
 323                 }
 324             }
 325
 326             // update reference markers
 327             refEntry.ID = al.RefID;
 328             lastRefID   = al.RefID;
 329             lastBin     = defaultValue;
 330         }
 331
 332         // if lastPosition greater than current alignment position - file not sorted properly
 333         else if ( lastPosition > al.Position ) {
 334             cerr << "BamStandardIndex ERROR: BAM file is not properly sorted by coordinate"
 335                  << ", aborting index creation"
 336                  << endl
 337                  << "At alignment: " << al.Name
 338                  << " : previous position " << lastPosition
 339                  << " > this alignment position " << al.Position
 340                  << " on reference id: " << al.RefID << endl;
 341             return false;
 342         }
 343
 344         // if alignment's ref ID is valid & its bin is not a 'leaf'
 345         if ( (al.RefID >= 0) && (al.Bin < 4681) )
 346             SaveLinearOffsetEntry(refEntry.LinearOffsets, al.Position, al.GetEndPosition(), lastOffset);
 347
 348         // changed to new BAI bin
 349         if ( al.Bin != lastBin ) {
 350
 351             // if not first bin on reference, save previous bin data
 352             if ( currentBin != defaultValue )
 353                 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
 354
 355             // update markers
 356             currentOffset = lastOffset;
 357             currentBin    = al.Bin;
 358             lastBin       = al.Bin;
 359             currentRefID  = al.RefID;
 360
 361             // if invalid RefID, break out
 362             if ( currentRefID < 0 )
 363                 break;
 364         }
 365
 366         // make sure that current file pointer is beyond lastOffset
 367         if ( m_reader->Tell() <= (int64_t)lastOffset ) {
 368             cerr << "BamStandardIndex ERROR: calculating offsets failed"
 369                  << ", aborting index creation" << endl;
 370             return false;
 371         }
 372
 373         // update lastOffset & lastPosition
 374         lastOffset   = m_reader->Tell();
 375         lastPosition = al.Position;
 376     }
 377
 378     // after finishing alignments, if any data was read, check:
 379     if ( currentRefID >= 0 ) {
 380
 381         // store last alignment chunk to its bin, then write last reference entry with data
 382         SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
 383         createdOk &= WriteReferenceEntry(refEntry);
 384
 385         // then write any empty references remaining at end of file
 386         for ( int i = currentRefID+1; i < numReferences; ++i ) {
 387             BaiReferenceEntry emptyEntry(i);
 388             createdOk &= WriteReferenceEntry(emptyEntry);
 389         }
 390     }
 391
 392     // rewind reader now that we're done building
 393     createdOk &= m_reader->Rewind();
 394
 395     // return result
 396     return createdOk;
 397 }
 398
 399 // returns format's file extension
 400 const string BamStandardIndex::Extension(void) {
 401     return BamStandardIndex::BAI_EXTENSION;
 402 }
 403
 404 bool BamStandardIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) {
 405
 406     // cannot calculate offsets if unknown/invalid reference ID requested
 407     if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() )
 408         return false;
 409
 410     // retrieve index summary for left bound reference
 411     const BaiReferenceSummary& refSummary = m_indexFileSummary.at(region.LeftRefID);
 412
 413     // set up region boundaries based on actual BamReader data
 414     uint32_t begin;
 415     uint32_t end;
 416     if ( !AdjustRegion(region, begin, end) ) {
 417         cerr << "BamStandardIndex ERROR: cannot calculate offsets on invalid region" << endl;
 418         return false;
 419     }
 420
 421     // retrieve all candidate bin IDs for region
 422     set<uint16_t> candidateBins;
 423     CalculateCandidateBins(begin, end, candidateBins);
 424
 425     // use reference's linear offsets to calculate the minimum offset
 426     // that must be considered to find overlap
 427     const uint64_t& minOffset = CalculateMinOffset(refSummary, begin);
 428
 429     // attempt to use reference summary, minOffset, & candidateBins to calculate offsets
 430     // no data should not be error
 431     vector<int64_t> offsets;
 432     if ( !CalculateCandidateOffsets(refSummary, minOffset, candidateBins, offsets) ) {
 433         cerr << "BamStandardIndex ERROR: could not calculate candidate offsets for requested region" << endl;
 434         return false;
 435     }
 436
 437     // if not candidate offsets are present in the indexed (most likely sparce coverage)
 438     // then silently bail
 439     if( offsets.size() == 0 ) {
 440         return false;
 441     }
 442
 443     // ensure that offsets are sorted before processing
 444     sort( offsets.begin(), offsets.end() );
 445
 446     // binary search for an overlapping block (may not be first one though)
 447     BamAlignment al;
 448     typedef vector<int64_t>::const_iterator OffsetConstIterator;
 449     OffsetConstIterator offsetFirst = offsets.begin();
 450     OffsetConstIterator offsetIter  = offsetFirst;
 451     OffsetConstIterator offsetLast  = offsets.end();
 452     iterator_traits<OffsetConstIterator>::difference_type count = distance(offsetFirst, offsetLast);
 453     iterator_traits<OffsetConstIterator>::difference_type step;
 454     while ( count > 0 ) {
 455         offsetIter = offsetFirst;
 456         step = count/2;
 457         advance(offsetIter, step);
 458
 459         // attempt seek to candidate offset
 460         const int64_t& candidateOffset = (*offsetIter);
 461         if ( !m_reader->Seek(candidateOffset) ) {
 462             cerr << "BamStandardIndex ERROR: could not jump"
 463                  << ", there was a problem seeking in BAM file" << endl;
 464             return false;
 465         }
 466
 467         // load first available alignment, setting flag to true if data exists
 468         *hasAlignmentsInRegion = m_reader->LoadNextAlignment(al);
 469
 470         // check alignment against region
 471         if ( al.GetEndPosition() < region.LeftPosition ) {
 472             offsetFirst = ++offsetIter;
 473             count -= step+1;
 474         } else count = step;
 475     }
 476
 477     // seek back to the offset before the 'current offset' (to cover overlaps)
 478     if ( offsetIter != offsets.begin() )
 479         --offsetIter;
 480     offset = (*offsetIter);
 481
 482     // return succes
 483     return true;
 484 }
 485
 486 // returns whether reference has alignments or no
 487 bool BamStandardIndex::HasAlignments(const int& referenceID) const {
 488     if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() )
 489         return false;
 490     const BaiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID);
 491     return ( refSummary.NumBins > 0 );
 492 }
 493
 494 bool BamStandardIndex::IsFileOpen(void) const {
 495     return ( m_indexStream != 0 );
 496 }
 497
 498 // attempts to use index data to jump to @region, returns success/fail
 499 // a "successful" jump indicates no error, but not whether this region has data
 500 //   * thus, the method sets a flag to indicate whether there are alignments
 501 //     available after the jump position
 502 bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
 503
 504     // clear out flag
 505     *hasAlignmentsInRegion = false;
 506
 507     // skip if reader is not valid or is not open
 508     if ( m_reader == 0 || !m_reader->IsOpen() )
 509         return false;
 510
 511     // calculate nearest offset to jump to
 512     int64_t offset;
 513     if ( !GetOffset(region, offset, hasAlignmentsInRegion) ) {
 514         cerr << "BamStandardIndex ERROR: could not jump"
 515              << ", unable to calculate offset for specified region" << endl;
 516         return false;
 517     }
 518
 519     // if region has alignments, return success/fail of seeking there
 520     if ( *hasAlignmentsInRegion )
 521         return m_reader->Seek(offset);
 522
 523     // otherwise, simply return true (but hasAlignmentsInRegion flag has been set to false)
 524     // (this is OK, BamReader will check this flag before trying to load data)
 525     return true;
 526 }
 527
 528 // loads existing data from file into memory
 529 bool BamStandardIndex::Load(const std::string& filename) {
 530
 531     // attempt open index file (read-only)
 532     if ( !OpenFile(filename, "rb") ) {
 533         cerr << "BamStandardIndex ERROR: could not open input index file: " << filename
 534              << ", aborting index load" << endl;
 535         return false;
 536     }
 537
 538     // if invalid format 'magic number', close & return failure
 539     if ( !CheckMagicNumber() ) {
 540         cerr << "BamStandardIndex ERROR: unexpected format for index file: " << filename
 541              << ", aborting index load" << endl;
 542         CloseFile();
 543         return false;
 544     }
 545
 546     // attempt to load index file summary, return success/failure
 547     if ( !SummarizeIndexFile() ) {
 548         cerr << "BamStandardIndex ERROR: could not generate a summary of index file " << filename
 549              << ", aborting index load" << endl;
 550         CloseFile();
 551         return false;
 552     }
 553
 554     // if we get here, index summary is loaded OK
 555     return true;
 556 }
 557
 558 uint64_t BamStandardIndex::LookupLinearOffset(const BaiReferenceSummary& refSummary, const int& index) {
 559
 560     // attempt seek to proper index file position
 561     const int64_t linearOffsetFilePosition = (int64_t)refSummary.FirstLinearOffsetFilePosition +
 562                                              index*BamStandardIndex::SIZEOF_LINEAROFFSET;
 563     if ( !Seek(linearOffsetFilePosition, SEEK_SET) )
 564         return 0;
 565
 566     // read linear offset from BAI file
 567     uint64_t linearOffset(0);
 568     if ( !ReadLinearOffset(linearOffset) )
 569         return 0;
 570     return linearOffset;
 571 }
 572
 573 void BamStandardIndex::MergeAlignmentChunks(BaiAlignmentChunkVector& chunks) {
 574
 575     // skip if chunks are empty, nothing to merge
 576     if ( chunks.empty() )
 577         return;
 578
 579     // set up merged alignment chunk container
 580     BaiAlignmentChunkVector mergedChunks;
 581     mergedChunks.push_back( chunks[0] );
 582
 583     // iterate over chunks
 584     int i = 0;
 585     BaiAlignmentChunkVector::iterator chunkIter = chunks.begin();
 586     BaiAlignmentChunkVector::iterator chunkEnd  = chunks.end();
 587     for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {
 588
 589         // get 'currentMergeChunk' based on numeric index
 590         BaiAlignmentChunk& currentMergeChunk = mergedChunks[i];
 591
 592         // get sourceChunk based on source vector iterator
 593         BaiAlignmentChunk& sourceChunk = (*chunkIter);
 594
 595         // if currentMergeChunk ends where sourceChunk starts, then merge the two
 596         if ( currentMergeChunk.Stop>>16 == sourceChunk.Start>>16 )
 597             currentMergeChunk.Stop = sourceChunk.Stop;
 598
 599         // otherwise
 600         else {
 601             // append sourceChunk after currentMergeChunk
 602             mergedChunks.push_back(sourceChunk);
 603
 604             // update i, so the next iteration will consider the
 605             // recently-appended sourceChunk as new mergeChunk candidate
 606             ++i;
 607         }
 608     }
 609
 610     // saved newly-merged chunks into (parameter) chunks
 611     chunks = mergedChunks;
 612 }
 613
 614 bool BamStandardIndex::OpenFile(const std::string& filename, const char* mode) {
 615
 616     // make sure any previous index file is closed
 617     CloseFile();
 618
 619     // attempt to open file
 620     m_indexStream = fopen(filename.c_str(), mode);
 621     return IsFileOpen();
 622 }
 623
 624 bool BamStandardIndex::ReadBinID(uint32_t& binId) {
 625     size_t elementsRead = 0;
 626     elementsRead += fread(&binId, sizeof(binId), 1, m_indexStream);
 627     if ( m_isBigEndian ) SwapEndian_32(binId);
 628     return ( elementsRead == 1 );
 629 }
 630
 631 bool BamStandardIndex::ReadBinIntoBuffer(uint32_t& binId, int32_t& numAlignmentChunks) {
 632
 633     bool readOk = true;
 634
 635     // read bin header
 636     readOk &= ReadBinID(binId);
 637     readOk &= ReadNumAlignmentChunks(numAlignmentChunks);
 638
 639     // read bin contents
 640     const unsigned int bytesRequested = numAlignmentChunks*BamStandardIndex::SIZEOF_ALIGNMENTCHUNK;
 641     readOk &= ReadIntoBuffer(bytesRequested);
 642
 643     // return success/failure
 644     return readOk;
 645 }
 646
 647 bool BamStandardIndex::ReadIntoBuffer(const unsigned int& bytesRequested) {
 648
 649     // ensure that our buffer is big enough for request
 650     BamStandardIndex::CheckBufferSize(m_buffer, m_bufferLength, bytesRequested);
 651
 652     // read from BAI file stream
 653     size_t bytesRead = fread( m_buffer, sizeof(char), bytesRequested, m_indexStream );
 654     return ( bytesRead == (size_t)bytesRequested );
 655 }
 656
 657 bool BamStandardIndex::ReadLinearOffset(uint64_t& linearOffset) {
 658     size_t elementsRead = 0;
 659     elementsRead += fread(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
 660     if ( m_isBigEndian ) SwapEndian_64(linearOffset);
 661     return ( elementsRead == 1 );
 662 }
 663
 664 bool BamStandardIndex::ReadNumAlignmentChunks(int& numAlignmentChunks) {
 665     size_t elementsRead = 0;
 666     elementsRead += fread(&numAlignmentChunks, sizeof(numAlignmentChunks), 1, m_indexStream);
 667     if ( m_isBigEndian ) SwapEndian_32(numAlignmentChunks);
 668     return ( elementsRead == 1 );
 669 }
 670
 671 bool BamStandardIndex::ReadNumBins(int& numBins) {
 672     size_t elementsRead = 0;
 673     elementsRead += fread(&numBins, sizeof(numBins), 1, m_indexStream);
 674     if ( m_isBigEndian ) SwapEndian_32(numBins);
 675     return ( elementsRead == 1 );
 676 }
 677
 678 bool BamStandardIndex::ReadNumLinearOffsets(int& numLinearOffsets) {
 679     size_t elementsRead = 0;
 680     elementsRead += fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, m_indexStream);
 681     if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets);
 682     return ( elementsRead == 1 );
 683 }
 684
 685 bool BamStandardIndex::ReadNumReferences(int& numReferences) {
 686     size_t elementsRead = 0;
 687     elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
 688     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 689     return ( elementsRead == 1 );
 690 }
 691
 692 void BamStandardIndex::ReserveForSummary(const int& numReferences) {
 693     m_indexFileSummary.clear();
 694     m_indexFileSummary.assign( numReferences, BaiReferenceSummary() );
 695 }
 696
 697 void BamStandardIndex::SaveAlignmentChunkToBin(BaiBinMap& binMap,
 698                                     const uint32_t& currentBin,
 699                                     const uint64_t& currentOffset,
 700                                     const uint64_t& lastOffset)
 701 {
 702     // create new alignment chunk
 703     BaiAlignmentChunk newChunk(currentOffset, lastOffset);
 704
 705
 706
 707     // if no entry exists yet for this bin, create one and store alignment chunk
 708     BaiBinMap::iterator binIter = binMap.find(currentBin);
 709     if ( binIter == binMap.end() ) {
 710         BaiAlignmentChunkVector newChunks;
 711         newChunks.push_back(newChunk);
 712         binMap.insert( pair<uint32_t, BaiAlignmentChunkVector>(currentBin, newChunks));
 713     }
 714
 715     // otherwise, just append alignment chunk
 716     else {
 717         BaiAlignmentChunkVector& binChunks = (*binIter).second;
 718         binChunks.push_back( newChunk );
 719     }
 720 }
 721
 722 void BamStandardIndex::SaveBinsSummary(const int& refId, const int& numBins) {
 723     BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
 724     refSummary.NumBins = numBins;
 725     refSummary.FirstBinFilePosition = Tell();
 726 }
 727
 728 void BamStandardIndex::SaveLinearOffsetEntry(BaiLinearOffsetVector& offsets,
 729                                              const int& alignmentStartPosition,
 730                                              const int& alignmentStopPosition,
 731                                              const uint64_t& lastOffset)
 732 {
 733     // get converted offsets
 734     const int beginOffset = alignmentStartPosition >> BamStandardIndex::BAM_LIDX_SHIFT;
 735     const int endOffset   = (alignmentStopPosition - 1) >> BamStandardIndex::BAM_LIDX_SHIFT;
 736
 737     // resize vector if necessary
 738     int oldSize = offsets.size();
 739     int newSize = endOffset + 1;
 740     if ( oldSize < newSize )
 741         offsets.resize(newSize, 0);
 742
 743     // store offset
 744     for( int i = beginOffset + 1; i <= endOffset; ++i ) {
 745         if ( offsets[i] == 0 )
 746             offsets[i] = lastOffset;
 747     }
 748 }
 749
 750 void BamStandardIndex::SaveLinearOffsetsSummary(const int& refId, const int& numLinearOffsets) {
 751     BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
 752     refSummary.NumLinearOffsets = numLinearOffsets;
 753     refSummary.FirstLinearOffsetFilePosition = Tell();
 754 }
 755
 756 // seek to position in index file stream
 757 bool BamStandardIndex::Seek(const int64_t& position, const int& origin) {
 758     return ( fseek64(m_indexStream, position, origin) == 0 );
 759 }
 760
 761 // change the index caching behavior
 762 void BamStandardIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) {
 763     m_cacheMode = mode;
 764     // do nothing else here ? cache mode will be ignored from now on, most likely
 765 }
 766
 767 bool BamStandardIndex::SkipBins(const int& numBins) {
 768     uint32_t binId;
 769     int32_t numAlignmentChunks;
 770     bool skippedOk = true;
 771     for (int i = 0; i < numBins; ++i)
 772         skippedOk &= ReadBinIntoBuffer(binId, numAlignmentChunks); // results & buffer ignored
 773     return skippedOk;
 774 }
 775
 776 bool BamStandardIndex::SkipLinearOffsets(const int& numLinearOffsets) {
 777     const unsigned int bytesRequested = numLinearOffsets*BamStandardIndex::SIZEOF_LINEAROFFSET;
 778     return ReadIntoBuffer(bytesRequested);
 779 }
 780
 781 void BamStandardIndex::SortLinearOffsets(BaiLinearOffsetVector& linearOffsets) {
 782     sort( linearOffsets.begin(), linearOffsets.end() );
 783 }
 784
 785 bool BamStandardIndex::SummarizeBins(BaiReferenceSummary& refSummary) {
 786
 787     // load number of bins
 788     int numBins;
 789     if ( !ReadNumBins(numBins) )
 790         return false;
 791
 792     // store bins summary for this reference
 793     refSummary.NumBins = numBins;
 794     refSummary.FirstBinFilePosition = Tell();
 795
 796     // attempt skip reference bins, return success/failure
 797     if ( !SkipBins(numBins) )
 798         return false;
 799
 800     // if we get here, bin summarized OK
 801     return true;
 802 }
 803
 804 bool BamStandardIndex::SummarizeIndexFile(void) {
 805
 806     // load number of reference sequences
 807     int numReferences;
 808     if ( !ReadNumReferences(numReferences) )
 809         return false;
 810
 811     // initialize file summary data
 812     ReserveForSummary(numReferences);
 813
 814     // iterate over reference entries
 815     bool loadedOk = true;
 816     BaiFileSummary::iterator summaryIter = m_indexFileSummary.begin();
 817     BaiFileSummary::iterator summaryEnd  = m_indexFileSummary.end();
 818     for ( int i = 0; summaryIter != summaryEnd; ++summaryIter, ++i )
 819         loadedOk &= SummarizeReference(*summaryIter);
 820
 821     // return result
 822     return loadedOk;
 823 }
 824
 825 bool BamStandardIndex::SummarizeLinearOffsets(BaiReferenceSummary& refSummary) {
 826
 827     // load number of linear offsets
 828     int numLinearOffsets;
 829     if ( !ReadNumLinearOffsets(numLinearOffsets) )
 830         return false;
 831
 832     // store bin summary data for this reference
 833     refSummary.NumLinearOffsets = numLinearOffsets;
 834     refSummary.FirstLinearOffsetFilePosition = Tell();
 835
 836     // skip linear offsets in index file
 837     if ( !SkipLinearOffsets(numLinearOffsets) )
 838         return false;
 839
 840     // if get here, linear offsets summarized OK
 841     return true;
 842 }
 843
 844 bool BamStandardIndex::SummarizeReference(BaiReferenceSummary& refSummary) {
 845
 846     bool loadedOk = true;
 847     loadedOk &= SummarizeBins(refSummary);
 848     loadedOk &= SummarizeLinearOffsets(refSummary);
 849     return loadedOk;
 850 }
 851
 852 // return position of file pointer in index file stream
 853 int64_t BamStandardIndex::Tell(void) const {
 854     return ftell64(m_indexStream);
 855 }
 856
 857 bool BamStandardIndex::WriteAlignmentChunk(const BaiAlignmentChunk& chunk) {
 858
 859     size_t elementsWritten = 0;
 860
 861     // localize alignment chunk offsets
 862     uint64_t start = chunk.Start;
 863     uint64_t stop  = chunk.Stop;
 864
 865     // swap endian-ness if necessary
 866     if ( m_isBigEndian ) {
 867         SwapEndian_64(start);
 868         SwapEndian_64(stop);
 869     }
 870
 871     // write to index file
 872     elementsWritten += fwrite(&start, sizeof(start), 1, m_indexStream);
 873     elementsWritten += fwrite(&stop,  sizeof(stop),  1, m_indexStream);
 874
 875     // return success/failure of write
 876     return ( elementsWritten == 2 );
 877 }
 878
 879 bool BamStandardIndex::WriteAlignmentChunks(BaiAlignmentChunkVector& chunks) {
 880
 881     // make sure chunks are merged (simplified) before writing & saving summary
 882     MergeAlignmentChunks(chunks);
 883
 884     size_t elementsWritten = 0;
 885
 886     // write chunks
 887     int32_t chunkCount = chunks.size();
 888     if ( m_isBigEndian ) SwapEndian_32(chunkCount);
 889     elementsWritten += fwrite(&chunkCount, sizeof(chunkCount), 1, m_indexStream);
 890
 891     // iterate over chunks
 892     bool chunksOk = true;
 893     BaiAlignmentChunkVector::const_iterator chunkIter = chunks.begin();
 894     BaiAlignmentChunkVector::const_iterator chunkEnd  = chunks.end();
 895     for ( ; chunkIter != chunkEnd; ++chunkIter )
 896         chunksOk &= WriteAlignmentChunk( (*chunkIter) );
 897
 898     // return success/failure of write
 899     return ( (elementsWritten == 1) && chunksOk );
 900 }
 901
 902 bool BamStandardIndex::WriteBin(const uint32_t& binId, BaiAlignmentChunkVector& chunks) {
 903
 904     size_t elementsWritten = 0;
 905
 906     // write BAM bin ID
 907     uint32_t binKey = binId;
 908     if ( m_isBigEndian ) SwapEndian_32(binKey);
 909     elementsWritten += fwrite(&binKey, sizeof(binKey), 1, m_indexStream);
 910
 911     // write bin's alignment chunks
 912     bool chunksOk = WriteAlignmentChunks(chunks);
 913
 914     // return success/failure of write
 915     return ( (elementsWritten == 1) && chunksOk );
 916 }
 917
 918 bool BamStandardIndex::WriteBins(const int& refId, BaiBinMap& bins) {
 919
 920     size_t elementsWritten = 0;
 921
 922     // write number of bins
 923     int32_t binCount = bins.size();
 924     if ( m_isBigEndian ) SwapEndian_32(binCount);
 925     elementsWritten += fwrite(&binCount, sizeof(binCount), 1, m_indexStream);
 926
 927     // save summary for reference's bins
 928     SaveBinsSummary(refId, bins.size());
 929
 930     // iterate over bins
 931     bool binsOk = true;
 932     BaiBinMap::iterator binIter = bins.begin();
 933     BaiBinMap::iterator binEnd  = bins.end();
 934     for ( ; binIter != binEnd; ++binIter )
 935         binsOk &= WriteBin( (*binIter).first, (*binIter).second );
 936
 937     // return success/failure of write
 938     return ( (elementsWritten == 1) && binsOk );
 939 }
 940
 941 bool BamStandardIndex::WriteHeader(void) {
 942
 943     size_t elementsWritten = 0;
 944
 945     // write magic number
 946     elementsWritten += fwrite(BamStandardIndex::BAI_MAGIC, sizeof(char), 4, m_indexStream);
 947
 948     // write number of reference sequences
 949     int32_t numReferences = m_indexFileSummary.size();
 950     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 951     elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, m_indexStream);
 952
 953     // return success/failure of write
 954     return (elementsWritten == 5);
 955 }
 956
 957 bool BamStandardIndex::WriteLinearOffsets(const int& refId, BaiLinearOffsetVector& linearOffsets) {
 958
 959     // make sure linear offsets are sorted before writing & saving summary
 960     SortLinearOffsets(linearOffsets);
 961
 962     size_t elementsWritten = 0;
 963
 964     // write number of linear offsets
 965     int32_t offsetCount = linearOffsets.size();
 966     if ( m_isBigEndian ) SwapEndian_32(offsetCount);
 967     elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, m_indexStream);
 968
 969     // save summary for reference's linear offsets
 970     SaveLinearOffsetsSummary(refId, linearOffsets.size());
 971
 972     // iterate over linear offsets
 973     BaiLinearOffsetVector::const_iterator offsetIter = linearOffsets.begin();
 974     BaiLinearOffsetVector::const_iterator offsetEnd  = linearOffsets.end();
 975     for ( ; offsetIter != offsetEnd; ++offsetIter ) {
 976
 977         // write linear offset
 978         uint64_t linearOffset = (*offsetIter);
 979         if ( m_isBigEndian ) SwapEndian_64(linearOffset);
 980         elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
 981     }
 982
 983     // return success/failure of write
 984     return ( elementsWritten == (size_t)(linearOffsets.size() + 1) );
 985 }
 986
 987 bool BamStandardIndex::WriteReferenceEntry(BaiReferenceEntry& refEntry) {
 988     bool refOk = true;
 989     refOk &= WriteBins(refEntry.ID, refEntry.Bins);
 990     refOk &= WriteLinearOffsets(refEntry.ID, refEntry.LinearOffsets);
 991     return refOk;
 992 }