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