src/api/internal/BamToolsIndex_p.cpp

   1 // ***************************************************************************
   2 // BamToolsIndex.cpp (c) 2010 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 6 October 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides index operations for the BamTools index format (".bti")
   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/BamToolsIndex_p.h>
  14 #include <api/internal/BgzfStream_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 <iostream>
  23 #include <iterator>
  24 #include <map>
  25 using namespace std;
  26
  27 // --------------------------------
  28 // static BamToolsIndex constants
  29 // --------------------------------
  30
  31 const uint32_t BamToolsIndex::DEFAULT_BLOCK_LENGTH = 1000;
  32 const string BamToolsIndex::BTI_EXTENSION     = ".bti";
  33 const char* const BamToolsIndex::BTI_MAGIC    = "BTI\1";
  34 const int BamToolsIndex::SIZEOF_BLOCK         = sizeof(int32_t)*2 + sizeof(int64_t);
  35
  36 // ----------------------------
  37 // RaiiWrapper implementation
  38 // ----------------------------
  39
  40 BamToolsIndex::RaiiWrapper::RaiiWrapper(void)
  41     : IndexStream(0)
  42 { }
  43
  44 BamToolsIndex::RaiiWrapper::~RaiiWrapper(void) {
  45     if ( IndexStream )
  46         fclose(IndexStream);
  47 }
  48
  49 // ------------------------------
  50 // BamToolsIndex implementation
  51 // ------------------------------
  52
  53 // ctor
  54 BamToolsIndex::BamToolsIndex(Internal::BamReaderPrivate* reader)
  55     : BamIndex(reader)
  56     , m_cacheMode(BamIndex::LimitedIndexCaching)
  57     , m_blockSize(BamToolsIndex::DEFAULT_BLOCK_LENGTH)
  58     , m_inputVersion(0)
  59     , m_outputVersion(BTI_1_2) // latest version - used for writing new index files
  60 {
  61     m_isBigEndian = BamTools::SystemIsBigEndian();
  62 }
  63
  64 // dtor
  65 BamToolsIndex::~BamToolsIndex(void) {
  66     CloseFile();
  67 }
  68
  69 void BamToolsIndex::CheckMagicNumber(void) {
  70
  71     // read magic number
  72     char magic[4];
  73     size_t elementsRead = fread(magic, sizeof(char), 4, Resources.IndexStream);
  74     if ( elementsRead != 4 )
  75         throw BamException("BamToolsIndex::CheckMagicNumber", "could not read BTI magic number");
  76
  77     // validate expected magic number
  78     if ( strncmp(magic, BamToolsIndex::BTI_MAGIC, 4) != 0 )
  79         throw BamException("BamToolsIndex::CheckMagicNumber", "invalid BTI magic number");
  80 }
  81
  82 // check index file version, return true if OK
  83 void BamToolsIndex::CheckVersion(void) {
  84
  85     // read version from file
  86     size_t elementsRead = fread(&m_inputVersion, sizeof(m_inputVersion), 1, Resources.IndexStream);
  87     if ( elementsRead != 1 )
  88         throw BamException("BamToolsIndex::CheckVersion", "could not read format version");
  89     if ( m_isBigEndian ) SwapEndian_32(m_inputVersion);
  90
  91     // if version is negative, or zero
  92     if ( m_inputVersion <= 0 )
  93         throw BamException("BamToolsIndex::CheckVersion", "invalid format version");
  94
  95     // if version is newer than can be supported by this version of bamtools
  96     else if ( m_inputVersion > m_outputVersion ) {
  97         const string message = "unsupported format: this index was created by a newer version of BamTools. "
  98                                "Update your local version of BamTools to use the index file.";
  99         throw BamException("BamToolsIndex::CheckVersion", message);
 100     }
 101
 102     // ------------------------------------------------------------------
 103     // check for deprecated, unsupported versions
 104     // (the format had to be modified to accomodate a particular bug fix)
 105
 106     else if ( (Version)m_inputVersion == BamToolsIndex::BTI_1_0 ) {
 107         const string message = "unsupported format: this version of the index may not properly handle "
 108                                "reference ends. Please run 'bamtools index -bti -in yourData.bam' to "
 109                                "generate an up-to-date, fixed BTI file.";
 110         throw BamException("BamToolsIndex::CheckVersion", message);
 111     }
 112
 113     else if ( (Version)m_inputVersion == BamToolsIndex::BTI_1_1 ) {
 114         const string message = "unsupported format: this version of the index may not properly handle "
 115                                "empty references. Please run 'bamtools index -bti -in yourData.bam to "
 116                                "generate an up-to-date, fixed BTI file.";
 117         throw BamException("BamToolsIndex::CheckVersion", message);
 118     }
 119 }
 120
 121 void BamToolsIndex::ClearReferenceEntry(BtiReferenceEntry& refEntry) {
 122     refEntry.ID = -1;
 123     refEntry.Blocks.clear();
 124 }
 125
 126 void BamToolsIndex::CloseFile(void) {
 127     if ( IsFileOpen() ) {
 128         fclose(Resources.IndexStream);
 129         Resources.IndexStream = 0;
 130     }
 131     m_indexFileSummary.clear();
 132 }
 133
 134 // builds index from associated BAM file & writes out to index file
 135 bool BamToolsIndex::Create(void) {
 136
 137     // skip if BamReader is invalid or not open
 138     if ( m_reader == 0 || !m_reader->IsOpen() ) {
 139         SetErrorString("BamToolsIndex::Create", "could not create index: reader is not open");
 140         return false;
 141     }
 142
 143     // rewind BamReader
 144     if ( !m_reader->Rewind() ) {
 145         const string readerError = m_reader->GetErrorString();
 146         const string message = "could not create index: \n\t" + readerError;
 147         SetErrorString("BamToolsIndex::Create", message);
 148         return false;
 149     }
 150
 151     try {
 152         // open new index file (read & write)
 153         string indexFilename = m_reader->Filename() + Extension();
 154         OpenFile(indexFilename, "w+b");
 155
 156         // initialize BtiFileSummary with number of references
 157         const int& numReferences = m_reader->GetReferenceCount();
 158         InitializeFileSummary(numReferences);
 159
 160         // intialize output file header
 161         WriteHeader();
 162
 163         // index building markers
 164         uint32_t currentBlockCount      = 0;
 165         int64_t currentAlignmentOffset  = m_reader->Tell();
 166         int32_t blockRefId              = -1;
 167         int32_t blockMaxEndPosition     = -1;
 168         int64_t blockStartOffset        = currentAlignmentOffset;
 169         int32_t blockStartPosition      = -1;
 170
 171         // plow through alignments, storing index entries
 172         BamAlignment al;
 173         BtiReferenceEntry refEntry;
 174         while ( m_reader->LoadNextAlignment(al) ) {
 175
 176             // if moved to new reference
 177             if ( al.RefID != blockRefId ) {
 178
 179                 // if first pass, check:
 180                 if ( currentBlockCount == 0 ) {
 181
 182                     // write any empty references up to (but not including) al.RefID
 183                     for ( int i = 0; i < al.RefID; ++i )
 184                         WriteReferenceEntry( BtiReferenceEntry(i) );
 185                 }
 186
 187                 // not first pass:
 188                 else {
 189
 190                     // store previous BTI block data in reference entry
 191                     BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 192                     refEntry.Blocks.push_back(block);
 193
 194                     // write reference entry, then clear
 195                     WriteReferenceEntry(refEntry);
 196                     ClearReferenceEntry(refEntry);
 197
 198                     // write any empty references between (but not including)
 199                     // the last blockRefID and current al.RefID
 200                     for ( int i = blockRefId+1; i < al.RefID; ++i )
 201                         WriteReferenceEntry( BtiReferenceEntry(i) );
 202
 203                     // reset block count
 204                     currentBlockCount = 0;
 205                 }
 206
 207                 // set ID for new reference entry
 208                 refEntry.ID = al.RefID;
 209             }
 210
 211             // if beginning of block, update counters
 212             if ( currentBlockCount == 0 ) {
 213                 blockRefId          = al.RefID;
 214                 blockStartOffset    = currentAlignmentOffset;
 215                 blockStartPosition  = al.Position;
 216                 blockMaxEndPosition = al.GetEndPosition();
 217             }
 218
 219             // increment block counter
 220             ++currentBlockCount;
 221
 222             // check end position
 223             int32_t alignmentEndPosition = al.GetEndPosition();
 224             if ( alignmentEndPosition > blockMaxEndPosition )
 225                 blockMaxEndPosition = alignmentEndPosition;
 226
 227             // if block is full, get offset for next block, reset currentBlockCount
 228             if ( currentBlockCount == m_blockSize ) {
 229
 230                 // store previous block data in reference entry
 231                 BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 232                 refEntry.Blocks.push_back(block);
 233
 234                 // update markers
 235                 blockStartOffset  = m_reader->Tell();
 236                 currentBlockCount = 0;
 237             }
 238
 239             // not the best name, but for the next iteration, this value will be the offset of the
 240             // *current* alignment. this is necessary because we won't know if this next alignment
 241             // is on a new reference until we actually read it
 242             currentAlignmentOffset = m_reader->Tell();
 243         }
 244
 245         // after finishing alignments, if any data was read, check:
 246         if ( blockRefId >= 0 ) {
 247
 248             // store last BTI block data in reference entry
 249             BtiBlock block(blockMaxEndPosition, blockStartOffset, blockStartPosition);
 250             refEntry.Blocks.push_back(block);
 251
 252             // write last reference entry, then clear
 253             WriteReferenceEntry(refEntry);
 254             ClearReferenceEntry(refEntry);
 255
 256             // then write any empty references remaining at end of file
 257             for ( int i = blockRefId+1; i < numReferences; ++i )
 258                 WriteReferenceEntry( BtiReferenceEntry(i) );
 259         }
 260
 261     } catch ( BamException& e ) {
 262         m_errorString = e.what();
 263         return false;
 264     }
 265
 266     // rewind BamReader
 267     if ( !m_reader->Rewind() ) {
 268         const string readerError = m_reader->GetErrorString();
 269         const string message = "could not create index: \n\t" + readerError;
 270         SetErrorString("BamToolsIndex::Create", message);
 271         return false;
 272     }
 273
 274     // return success
 275     return true;
 276 }
 277
 278 // returns format's file extension
 279 const std::string BamToolsIndex::Extension(void) {
 280     return BamToolsIndex::BTI_EXTENSION;
 281 }
 282
 283 void BamToolsIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) {
 284
 285     // return false ref ID is not a valid index in file summary data
 286     if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() )
 287         throw BamException("BamToolsIndex::GetOffset", "invalid region requested");
 288
 289     // retrieve reference index data for left bound reference
 290     BtiReferenceEntry refEntry(region.LeftRefID);
 291     ReadReferenceEntry(refEntry);
 292
 293     // binary search for an overlapping block (may not be first one though)
 294     bool found = false;
 295     typedef BtiBlockVector::const_iterator BtiBlockConstIterator;
 296     BtiBlockConstIterator blockFirst = refEntry.Blocks.begin();
 297     BtiBlockConstIterator blockIter  = blockFirst;
 298     BtiBlockConstIterator blockLast  = refEntry.Blocks.end();
 299     iterator_traits<BtiBlockConstIterator>::difference_type count = distance(blockFirst, blockLast);
 300     iterator_traits<BtiBlockConstIterator>::difference_type step;
 301     while ( count > 0 ) {
 302         blockIter = blockFirst;
 303         step = count/2;
 304         advance(blockIter, step);
 305
 306         const BtiBlock& block = (*blockIter);
 307         if ( block.StartPosition <= region.RightPosition ) {
 308             if ( block.MaxEndPosition >= region.LeftPosition ) {
 309                 offset = block.StartOffset;
 310                 break;
 311             }
 312             blockFirst = ++blockIter;
 313             count -= step+1;
 314         }
 315         else count = step;
 316     }
 317
 318     // if we didn't search "off the end" of the blocks
 319     if ( blockIter != blockLast ) {
 320
 321         // "walk back" until we've gone too far
 322         while ( blockIter != blockFirst ) {
 323             const BtiBlock& currentBlock = (*blockIter);
 324
 325             --blockIter;
 326             const BtiBlock& previousBlock = (*blockIter);
 327             if ( previousBlock.MaxEndPosition < region.LeftPosition ) {
 328                 offset = currentBlock.StartOffset;
 329                 found = true;
 330                 break;
 331             }
 332         }
 333
 334         // if we walked all the way to first block, just return that and let the reader's
 335         // region overlap parsing do the rest
 336         if ( blockIter == blockFirst ) {
 337             const BtiBlock& block = (*blockIter);
 338             offset = block.StartOffset;
 339             found = true;
 340         }
 341     }
 342
 343
 344     // sets to false if blocks container is empty, or if no matching block could be found
 345     *hasAlignmentsInRegion = found;
 346 }
 347
 348 // returns whether reference has alignments or no
 349 bool BamToolsIndex::HasAlignments(const int& referenceID) const {
 350     if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() )
 351         return false;
 352     const BtiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID);
 353     return ( refSummary.NumBlocks > 0 );
 354 }
 355
 356 // pre-allocates space for each reference's summary data
 357 void BamToolsIndex::InitializeFileSummary(const int& numReferences) {
 358     m_indexFileSummary.clear();
 359     for ( int i = 0; i < numReferences; ++i )
 360         m_indexFileSummary.push_back( BtiReferenceSummary() );
 361 }
 362
 363 // returns true if the index stream is open
 364 bool BamToolsIndex::IsFileOpen(void) const {
 365     return ( Resources.IndexStream != 0 );
 366 }
 367
 368 // attempts to use index data to jump to @region, returns success/fail
 369 // a "successful" jump indicates no error, but not whether this region has data
 370 //   * thus, the method sets a flag to indicate whether there are alignments
 371 //     available after the jump position
 372 bool BamToolsIndex::Jump(const BamTools::BamRegion& region, bool* hasAlignmentsInRegion) {
 373
 374     // clear flag
 375     *hasAlignmentsInRegion = false;
 376
 377     // skip if invalid reader or not open
 378     if ( m_reader == 0 || !m_reader->IsOpen() ) {
 379         SetErrorString("BamToolsIndex::Jump", "could not jump: reader is not open");
 380         return false;
 381     }
 382
 383     // make sure left-bound position is valid
 384     const RefVector& references = m_reader->GetReferenceData();
 385     if ( region.LeftPosition > references.at(region.LeftRefID).RefLength ) {
 386         SetErrorString("BamToolsIndex::Jump", "could not create index: invalid region requested");
 387         return false;
 388     }
 389
 390     // calculate nearest offset to jump to
 391     int64_t offset;
 392     try {
 393         GetOffset(region, offset, hasAlignmentsInRegion);
 394     } catch ( BamException& e ) {
 395         m_errorString = e.what();
 396         return false;
 397     }
 398
 399     // return success/failure of seek
 400     return m_reader->Seek(offset);
 401 }
 402
 403 // loads existing data from file into memory
 404 bool BamToolsIndex::Load(const std::string& filename) {
 405
 406     try {
 407
 408         // attempt to open file (read-only)
 409         OpenFile(filename, "rb");
 410
 411         // load metadata & generate in-memory summary
 412         LoadHeader();
 413         LoadFileSummary();
 414
 415         // return success
 416         return true;
 417
 418     } catch ( BamException& e ) {
 419         m_errorString = e.what();
 420         return false;
 421     }
 422 }
 423
 424 void BamToolsIndex::LoadFileSummary(void) {
 425
 426     // load number of reference sequences
 427     int numReferences;
 428     LoadNumReferences(numReferences);
 429
 430     // initialize file summary data
 431     InitializeFileSummary(numReferences);
 432
 433     // load summary for each reference
 434     BtiFileSummary::iterator summaryIter = m_indexFileSummary.begin();
 435     BtiFileSummary::iterator summaryEnd  = m_indexFileSummary.end();
 436     for ( ; summaryIter != summaryEnd; ++summaryIter )
 437         LoadReferenceSummary(*summaryIter);
 438 }
 439
 440 void BamToolsIndex::LoadHeader(void) {
 441
 442     // check BTI file metadata
 443     CheckMagicNumber();
 444     CheckVersion();
 445
 446     // use file's BTI block size to set member variable
 447     const size_t elementsRead = fread(&m_blockSize, sizeof(m_blockSize), 1, Resources.IndexStream);
 448     if ( m_isBigEndian ) SwapEndian_32(m_blockSize);
 449     if ( elementsRead != 1 )
 450         throw BamException("BamToolsIndex::LoadHeader", "could not read BTI block size");
 451 }
 452
 453 void BamToolsIndex::LoadNumBlocks(int& numBlocks) {
 454     const size_t elementsRead = fread(&numBlocks, sizeof(numBlocks), 1, Resources.IndexStream);
 455     if ( m_isBigEndian ) SwapEndian_32(numBlocks);
 456     if ( elementsRead != 1 )
 457         throw BamException("BamToolsIndex::LoadNumBlocks", "could not read number of BTI blocks");
 458 }
 459
 460 void BamToolsIndex::LoadNumReferences(int& numReferences) {
 461     const size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, Resources.IndexStream);
 462     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 463     if ( elementsRead != 1 )
 464         throw BamException("BamToolsIndex::LoadNumReferences", "could not read number of references");
 465 }
 466
 467 void BamToolsIndex::LoadReferenceSummary(BtiReferenceSummary& refSummary) {
 468
 469     // load number of blocks
 470     int numBlocks;
 471     LoadNumBlocks(numBlocks);
 472
 473     // store block summary data for this reference
 474     refSummary.NumBlocks = numBlocks;
 475     refSummary.FirstBlockFilePosition = Tell();
 476
 477     // skip reference's blocks
 478     SkipBlocks(numBlocks);
 479 }
 480
 481 void BamToolsIndex::OpenFile(const std::string& filename, const char* mode) {
 482
 483     // make sure any previous index file is closed
 484     CloseFile();
 485
 486     // attempt to open file
 487     Resources.IndexStream = fopen(filename.c_str(), mode);
 488     if ( !IsFileOpen() ) {
 489         const string message = string("could not open file: ") + filename;
 490         throw BamException("BamToolsIndex::OpenFile", message);
 491     }
 492 }
 493
 494 void BamToolsIndex::ReadBlock(BtiBlock& block) {
 495
 496     // read in block data members
 497     size_t elementsRead = 0;
 498     elementsRead += fread(&block.MaxEndPosition, sizeof(block.MaxEndPosition), 1, Resources.IndexStream);
 499     elementsRead += fread(&block.StartOffset,    sizeof(block.StartOffset),    1, Resources.IndexStream);
 500     elementsRead += fread(&block.StartPosition,  sizeof(block.StartPosition),  1, Resources.IndexStream);
 501
 502     // swap endian-ness if necessary
 503     if ( m_isBigEndian ) {
 504         SwapEndian_32(block.MaxEndPosition);
 505         SwapEndian_64(block.StartOffset);
 506         SwapEndian_32(block.StartPosition);
 507     }
 508
 509     if ( elementsRead != 3 )
 510         throw BamException("BamToolsIndex::ReadBlock", "could not read block");
 511 }
 512
 513 void BamToolsIndex::ReadBlocks(const BtiReferenceSummary& refSummary, BtiBlockVector& blocks) {
 514
 515     // prep blocks container
 516     blocks.clear();
 517     blocks.reserve(refSummary.NumBlocks);
 518
 519     // skip to first block entry
 520     Seek( refSummary.FirstBlockFilePosition, SEEK_SET );
 521
 522     // read & store block entries
 523     BtiBlock block;
 524     for ( int i = 0; i < refSummary.NumBlocks; ++i ) {
 525         ReadBlock(block);
 526         blocks.push_back(block);
 527     }
 528 }
 529
 530 void BamToolsIndex::ReadReferenceEntry(BtiReferenceEntry& refEntry) {
 531
 532     // return false if refId not valid index in file summary structure
 533     if ( refEntry.ID < 0 || refEntry.ID >= (int)m_indexFileSummary.size() )
 534         throw BamException("BamToolsIndex::ReadReferenceEntry", "invalid reference requested");
 535
 536     // use index summary to assist reading the reference's BTI blocks
 537     const BtiReferenceSummary& refSummary = m_indexFileSummary.at(refEntry.ID);
 538     ReadBlocks(refSummary, refEntry.Blocks);
 539 }
 540
 541 void BamToolsIndex::Seek(const int64_t& position, const int& origin) {
 542     if ( fseek64(Resources.IndexStream, position, origin) != 0 )
 543         throw BamException("BamToolsIndex::Seek", "could not seek in BAI file");
 544 }
 545
 546 // change the index caching behavior
 547 void BamToolsIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) {
 548     m_cacheMode = mode;
 549     // do nothing else here ? cache mode will be ignored from now on, most likely
 550 }
 551
 552 void BamToolsIndex::SkipBlocks(const int& numBlocks) {
 553     Seek( numBlocks*BamToolsIndex::SIZEOF_BLOCK, SEEK_CUR );
 554 }
 555
 556 int64_t BamToolsIndex::Tell(void) const {
 557     return ftell64(Resources.IndexStream);
 558 }
 559
 560 void BamToolsIndex::WriteBlock(const BtiBlock& block) {
 561
 562     // copy entry data
 563     int32_t maxEndPosition = block.MaxEndPosition;
 564     int64_t startOffset    = block.StartOffset;
 565     int32_t startPosition  = block.StartPosition;
 566
 567     // swap endian-ness if necessary
 568     if ( m_isBigEndian ) {
 569         SwapEndian_32(maxEndPosition);
 570         SwapEndian_64(startOffset);
 571         SwapEndian_32(startPosition);
 572     }
 573
 574     // write the reference index entry
 575     size_t elementsWritten = 0;
 576     elementsWritten += fwrite(&maxEndPosition, sizeof(maxEndPosition), 1, Resources.IndexStream);
 577     elementsWritten += fwrite(&startOffset,    sizeof(startOffset),    1, Resources.IndexStream);
 578     elementsWritten += fwrite(&startPosition,  sizeof(startPosition),  1, Resources.IndexStream);
 579     if ( elementsWritten != 3 )
 580         throw BamException("BamToolsIndex::WriteBlock", "could not write BTI block");
 581 }
 582
 583 void BamToolsIndex::WriteBlocks(const BtiBlockVector& blocks) {
 584     BtiBlockVector::const_iterator blockIter = blocks.begin();
 585     BtiBlockVector::const_iterator blockEnd  = blocks.end();
 586     for ( ; blockIter != blockEnd; ++blockIter )
 587         WriteBlock(*blockIter);
 588 }
 589
 590 void BamToolsIndex::WriteHeader(void) {
 591
 592     size_t elementsWritten = 0;
 593
 594     // write BTI index format 'magic number'
 595     elementsWritten += fwrite(BamToolsIndex::BTI_MAGIC, 1, 4, Resources.IndexStream);
 596
 597     // write BTI index format version
 598     int32_t currentVersion = (int32_t)m_outputVersion;
 599     if ( m_isBigEndian ) SwapEndian_32(currentVersion);
 600     elementsWritten += fwrite(&currentVersion, sizeof(currentVersion), 1, Resources.IndexStream);
 601
 602     // write block size
 603     uint32_t blockSize = m_blockSize;
 604     if ( m_isBigEndian ) SwapEndian_32(blockSize);
 605     elementsWritten += fwrite(&blockSize, sizeof(blockSize), 1, Resources.IndexStream);
 606
 607     // write number of references
 608     int32_t numReferences = m_indexFileSummary.size();
 609     if ( m_isBigEndian ) SwapEndian_32(numReferences);
 610     elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, Resources.IndexStream);
 611
 612     if ( elementsWritten != 7 )
 613         throw BamException("BamToolsIndex::WriteHeader", "could not write BTI header");
 614 }
 615
 616 void BamToolsIndex::WriteReferenceEntry(const BtiReferenceEntry& refEntry) {
 617
 618     // write number of blocks this reference
 619     uint32_t numBlocks = refEntry.Blocks.size();
 620     if ( m_isBigEndian ) SwapEndian_32(numBlocks);
 621     const size_t elementsWritten = fwrite(&numBlocks, sizeof(numBlocks), 1, Resources.IndexStream);
 622     if ( elementsWritten != 1 )
 623         throw BamException("BamToolsIndex::WriteReferenceEntry", "could not write number of blocks");
 624
 625     // write actual block entries
 626     WriteBlocks(refEntry.Blocks);
 627 }