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