src/api/internal/index/BamStandardIndex_p.cpp

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