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