BamReader.cpp

   1 // ***************************************************************************\r
   2 // BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg\r
   3 // Marth Lab, Department of Biology, Boston College\r
   4 // All rights reserved.\r
   5 // ---------------------------------------------------------------------------\r
   6 // Last modified: 29 March 2010 (DB)\r
   7 // ---------------------------------------------------------------------------\r
   8 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
   9 // Institute.\r
  10 // ---------------------------------------------------------------------------\r
  11 // Provides the basic functionality for reading BAM files\r
  12 // ***************************************************************************\r
  13 \r
  14 // C++ includes\r
  15 #include <algorithm>\r
  16 #include <iterator>\r
  17 #include <string>\r
  18 #include <vector>\r
  19 \r
  20 // BamTools includes\r
  21 #include "BGZF.h"\r
  22 #include "BamReader.h"\r
  23 using namespace BamTools;\r
  24 using namespace std;\r
  25 \r
  26 struct BamReader::BamReaderPrivate {\r
  27 \r
  28     // -------------------------------\r
  29     // data members\r
  30     // -------------------------------\r
  31 \r
  32     // general data\r
  33     BgzfData  mBGZF;\r
  34     string    HeaderText;\r
  35     BamIndex  Index;\r
  36     RefVector References;\r
  37     bool      IsIndexLoaded;\r
  38     int64_t   AlignmentsBeginOffset;\r
  39     string    Filename;\r
  40     string    IndexFilename;\r
  41     \r
  42     bool IsBigEndian;\r
  43 \r
  44     // user-specified region values\r
  45     bool IsRegionSpecified;\r
  46     int  CurrentRefID;\r
  47     int  CurrentLeft;\r
  48 \r
  49     // BAM character constants\r
  50     const char* DNA_LOOKUP;\r
  51     const char* CIGAR_LOOKUP;\r
  52 \r
  53     // -------------------------------\r
  54     // constructor & destructor\r
  55     // -------------------------------\r
  56     BamReaderPrivate(void);\r
  57     ~BamReaderPrivate(void);\r
  58 \r
  59     // -------------------------------\r
  60     // "public" interface\r
  61     // -------------------------------\r
  62 \r
  63     // flie operations\r
  64     void Close(void);\r
  65     bool Jump(int refID, int position = 0);\r
  66     void Open(const string& filename, const string& indexFilename = "");\r
  67     bool Rewind(void);\r
  68 \r
  69     // access alignment data\r
  70     bool GetNextAlignment(BamAlignment& bAlignment);\r
  71 \r
  72     // access auxiliary data\r
  73     const string GetHeaderText(void) const;\r
  74     const int GetReferenceCount(void) const;\r
  75     const RefVector GetReferenceData(void) const;\r
  76     const int GetReferenceID(const string& refName) const;\r
  77 \r
  78     // index operations\r
  79     bool CreateIndex(void);\r
  80 \r
  81     // -------------------------------\r
  82     // internal methods\r
  83     // -------------------------------\r
  84 \r
  85     // *** reading alignments and auxiliary data *** //\r
  86 \r
  87     // calculate bins that overlap region ( left to reference end for now )\r
  88     int BinsFromRegion(int refID, int left, uint16_t[MAX_BIN]);\r
  89     // calculates alignment end position based on starting position and provided CIGAR operations\r
  90     int CalculateAlignmentEnd(const int& position, const std::vector<CigarOp>& cigarData);\r
  91     // calculate file offset for first alignment chunk overlapping 'left'\r
  92     int64_t GetOffset(int refID, int left);\r
  93     // checks to see if alignment overlaps current region\r
  94     bool IsOverlap(BamAlignment& bAlignment);\r
  95     // retrieves header text from BAM file\r
  96     void LoadHeaderData(void);\r
  97     // retrieves BAM alignment under file pointer\r
  98     bool LoadNextAlignment(BamAlignment& bAlignment);\r
  99     // builds reference data structure from BAM file\r
 100     void LoadReferenceData(void);\r
 101 \r
 102     // *** index file handling *** //\r
 103 \r
 104     // calculates index for BAM file\r
 105     bool BuildIndex(void);\r
 106     // clear out inernal index data structure\r
 107     void ClearIndex(void);\r
 108     // saves BAM bin entry for index\r
 109     void InsertBinEntry(BamBinMap& binMap, const uint32_t& saveBin, const uint64_t& saveOffset, const uint64_t& lastOffset);\r
 110     // saves linear offset entry for index\r
 111     void InsertLinearOffset(LinearOffsetVector& offsets, const BamAlignment& bAlignment, const uint64_t& lastOffset);\r
 112     // loads index from BAM index file\r
 113     bool LoadIndex(void);\r
 114     // simplifies index by merging 'chunks'\r
 115     void MergeChunks(void);\r
 116     // round-up 32-bit integer to next power-of-2\r
 117     void Roundup32(int& value);\r
 118     // saves index to BAM index file\r
 119     bool WriteIndex(void);\r
 120 };\r
 121 \r
 122 // -----------------------------------------------------\r
 123 // BamReader implementation (wrapper around BRPrivate)\r
 124 // -----------------------------------------------------\r
 125 \r
 126 // constructor\r
 127 BamReader::BamReader(void) {\r
 128     d = new BamReaderPrivate;\r
 129 }\r
 130 \r
 131 // destructor\r
 132 BamReader::~BamReader(void) {\r
 133     delete d;\r
 134     d = 0;\r
 135 }\r
 136 \r
 137 // file operations\r
 138 void BamReader::Close(void) { d->Close(); }\r
 139 bool BamReader::Jump(int refID, int position) { return d->Jump(refID, position); }\r
 140 void BamReader::Open(const string& filename, const string& indexFilename) { d->Open(filename, indexFilename); }\r
 141 bool BamReader::Rewind(void) { return d->Rewind(); }\r
 142 \r
 143 // access alignment data\r
 144 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }\r
 145 \r
 146 // access auxiliary data\r
 147 const string    BamReader::GetHeaderText(void) const { return d->HeaderText; }\r
 148 const int       BamReader::GetReferenceCount(void) const { return d->References.size(); }\r
 149 const RefVector BamReader::GetReferenceData(void) const { return d->References; }\r
 150 const int       BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }\r
 151 \r
 152 // index operations\r
 153 bool BamReader::CreateIndex(void) { return d->CreateIndex(); }\r
 154 \r
 155 // -----------------------------------------------------\r
 156 // BamReaderPrivate implementation\r
 157 // -----------------------------------------------------\r
 158 \r
 159 // constructor\r
 160 BamReader::BamReaderPrivate::BamReaderPrivate(void)\r
 161     : IsIndexLoaded(false)\r
 162     , AlignmentsBeginOffset(0)\r
 163     , IsRegionSpecified(false)\r
 164     , CurrentRefID(0)\r
 165     , CurrentLeft(0)\r
 166     , DNA_LOOKUP("=ACMGRSVTWYHKDBN")\r
 167     , CIGAR_LOOKUP("MIDNSHP")\r
 168 { \r
 169     IsBigEndian = SystemIsBigEndian();\r
 170 }\r
 171 \r
 172 // destructor\r
 173 BamReader::BamReaderPrivate::~BamReaderPrivate(void) {\r
 174     Close();\r
 175 }\r
 176 \r
 177 // calculate bins that overlap region ( left to reference end for now )\r
 178 int BamReader::BamReaderPrivate::BinsFromRegion(int refID, int left, uint16_t list[MAX_BIN]) {\r
 179 \r
 180     // get region boundaries\r
 181     uint32_t begin = (unsigned int)left;\r
 182     uint32_t end   = (unsigned int)References.at(refID).RefLength - 1;\r
 183 \r
 184     // initialize list, bin '0' always a valid bin\r
 185     int i = 0;\r
 186     list[i++] = 0;\r
 187 \r
 188     // get rest of bins that contain this region\r
 189     unsigned int k;\r
 190     for (k =    1 + (begin>>26); k <=    1 + (end>>26); ++k) { list[i++] = k; }\r
 191     for (k =    9 + (begin>>23); k <=    9 + (end>>23); ++k) { list[i++] = k; }\r
 192     for (k =   73 + (begin>>20); k <=   73 + (end>>20); ++k) { list[i++] = k; }\r
 193     for (k =  585 + (begin>>17); k <=  585 + (end>>17); ++k) { list[i++] = k; }\r
 194     for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { list[i++] = k; }\r
 195 \r
 196     // return number of bins stored\r
 197     return i;\r
 198 }\r
 199 \r
 200 // populates BAM index data structure from BAM file data\r
 201 bool BamReader::BamReaderPrivate::BuildIndex(void) {\r
 202 \r
 203     // check to be sure file is open\r
 204     if (!mBGZF.IsOpen) { return false; }\r
 205 \r
 206     // move file pointer to beginning of alignments\r
 207     Rewind();\r
 208 \r
 209     // get reference count, reserve index space\r
 210     int numReferences = References.size();\r
 211     for ( int i = 0; i < numReferences; ++i ) {\r
 212         Index.push_back(ReferenceIndex());\r
 213     }\r
 214 \r
 215     // sets default constant for bin, ID, offset, coordinate variables\r
 216     const uint32_t defaultValue = 0xffffffffu;\r
 217 \r
 218     // bin data\r
 219     uint32_t saveBin(defaultValue);\r
 220     uint32_t lastBin(defaultValue);\r
 221 \r
 222     // reference ID data\r
 223     int32_t saveRefID(defaultValue);\r
 224     int32_t lastRefID(defaultValue);\r
 225 \r
 226     // offset data\r
 227     uint64_t saveOffset = mBGZF.Tell();\r
 228     uint64_t lastOffset = saveOffset;\r
 229 \r
 230     // coordinate data\r
 231     int32_t lastCoordinate = defaultValue;\r
 232 \r
 233     BamAlignment bAlignment;\r
 234     while( GetNextAlignment(bAlignment) ) {\r
 235 \r
 236         // change of chromosome, save ID, reset bin\r
 237         if ( lastRefID != bAlignment.RefID ) {\r
 238             lastRefID = bAlignment.RefID;\r
 239             lastBin   = defaultValue;\r
 240         }\r
 241 \r
 242         // if lastCoordinate greater than BAM position - file not sorted properly\r
 243         else if ( lastCoordinate > bAlignment.Position ) {\r
 244             printf("BAM file not properly sorted:\n");\r
 245             printf("Alignment %s : %d > %d on reference (id = %d)", bAlignment.Name.c_str(), lastCoordinate, bAlignment.Position, bAlignment.RefID);\r
 246             exit(1);\r
 247         }\r
 248 \r
 249         // if valid reference && BAM bin spans some minimum cutoff (smaller bin ids span larger regions)\r
 250         if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) {\r
 251 \r
 252             // save linear offset entry (matched to BAM entry refID)\r
 253             ReferenceIndex& refIndex = Index.at(bAlignment.RefID);\r
 254             LinearOffsetVector& offsets = refIndex.Offsets;\r
 255             InsertLinearOffset(offsets, bAlignment, lastOffset);\r
 256         }\r
 257 \r
 258         // if current BamAlignment bin != lastBin, "then possibly write the binning index"\r
 259         if ( bAlignment.Bin != lastBin ) {\r
 260 \r
 261             // if not first time through\r
 262             if ( saveBin != defaultValue ) {\r
 263 \r
 264                 // save Bam bin entry\r
 265                 ReferenceIndex& refIndex = Index.at(saveRefID);\r
 266                 BamBinMap& binMap = refIndex.Bins;\r
 267                 InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);\r
 268             }\r
 269 \r
 270             // update saveOffset\r
 271             saveOffset = lastOffset;\r
 272 \r
 273             // update bin values\r
 274             saveBin = bAlignment.Bin;\r
 275             lastBin = bAlignment.Bin;\r
 276 \r
 277             // update saveRefID\r
 278             saveRefID = bAlignment.RefID;\r
 279 \r
 280             // if invalid RefID, break out (why?)\r
 281             if ( saveRefID < 0 ) { break; }\r
 282         }\r
 283 \r
 284         // make sure that current file pointer is beyond lastOffset\r
 285         if ( mBGZF.Tell() <= (int64_t)lastOffset  ) {\r
 286             printf("Error in BGZF offsets.\n");\r
 287             exit(1);\r
 288         }\r
 289 \r
 290         // update lastOffset\r
 291         lastOffset = mBGZF.Tell();\r
 292 \r
 293         // update lastCoordinate\r
 294         lastCoordinate = bAlignment.Position;\r
 295     }\r
 296 \r
 297     // save any leftover BAM data (as long as refID is valid)\r
 298     if ( saveRefID >= 0 ) {\r
 299         // save Bam bin entry\r
 300         ReferenceIndex& refIndex = Index.at(saveRefID);\r
 301         BamBinMap& binMap = refIndex.Bins;\r
 302         InsertBinEntry(binMap, saveBin, saveOffset, lastOffset);\r
 303     }\r
 304 \r
 305     // simplify index by merging chunks\r
 306     MergeChunks();\r
 307 \r
 308     // iterate over references\r
 309     BamIndex::iterator indexIter = Index.begin();\r
 310     BamIndex::iterator indexEnd  = Index.end();\r
 311     for ( int i = 0; indexIter != indexEnd; ++indexIter, ++i ) {\r
 312 \r
 313         // get reference index data\r
 314         ReferenceIndex& refIndex = (*indexIter);\r
 315         BamBinMap& binMap = refIndex.Bins;\r
 316         LinearOffsetVector& offsets = refIndex.Offsets;\r
 317 \r
 318         // store whether reference has alignments or no\r
 319         References[i].RefHasAlignments = ( binMap.size() > 0 );\r
 320 \r
 321         // sort linear offsets\r
 322         sort(offsets.begin(), offsets.end());\r
 323     }\r
 324 \r
 325 \r
 326     // rewind file pointer to beginning of alignments, return success/fail\r
 327     return Rewind();\r
 328 }\r
 329 \r
 330 // calculates alignment end position based on starting position and provided CIGAR operations\r
 331 int BamReader::BamReaderPrivate::CalculateAlignmentEnd(const int& position, const vector<CigarOp>& cigarData) {\r
 332 \r
 333     // initialize alignment end to starting position\r
 334     int alignEnd = position;\r
 335 \r
 336     // iterate over cigar operations\r
 337     vector<CigarOp>::const_iterator cigarIter = cigarData.begin();\r
 338     vector<CigarOp>::const_iterator cigarEnd  = cigarData.end();\r
 339     for ( ; cigarIter != cigarEnd; ++cigarIter) {\r
 340         char cigarType = (*cigarIter).Type;\r
 341         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {\r
 342             alignEnd += (*cigarIter).Length;\r
 343         }\r
 344     }\r
 345     return alignEnd;\r
 346 }\r
 347 \r
 348 \r
 349 // clear index data structure\r
 350 void BamReader::BamReaderPrivate::ClearIndex(void) {\r
 351     Index.clear(); // sufficient ??\r
 352 }\r
 353 \r
 354 // closes the BAM file\r
 355 void BamReader::BamReaderPrivate::Close(void) {\r
 356     mBGZF.Close();\r
 357     ClearIndex();\r
 358     HeaderText.clear();\r
 359     IsRegionSpecified = false;\r
 360 }\r
 361 \r
 362 // create BAM index from BAM file (keep structure in memory) and write to default index output file\r
 363 bool BamReader::BamReaderPrivate::CreateIndex(void) {\r
 364 \r
 365     // clear out index\r
 366     ClearIndex();\r
 367 \r
 368     // build (& save) index from BAM file\r
 369     bool ok = true;\r
 370     ok &= BuildIndex();\r
 371     ok &= WriteIndex();\r
 372 \r
 373     // return success/fail\r
 374     return ok;\r
 375 }\r
 376 \r
 377 // returns RefID for given RefName (returns References.size() if not found)\r
 378 const int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {\r
 379 \r
 380     // retrieve names from reference data\r
 381     vector<string> refNames;\r
 382     RefVector::const_iterator refIter = References.begin();\r
 383     RefVector::const_iterator refEnd  = References.end();\r
 384     for ( ; refIter != refEnd; ++refIter) {\r
 385         refNames.push_back( (*refIter).RefName );\r
 386     }\r
 387 \r
 388     // return 'index-of' refName ( if not found, returns refNames.size() )\r
 389     return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));\r
 390 }\r
 391 \r
 392 // get next alignment (from specified region, if given)\r
 393 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {\r
 394 \r
 395     // if valid alignment available\r
 396     if ( LoadNextAlignment(bAlignment) ) {\r
 397 \r
 398         // if region not specified, return success\r
 399         if ( !IsRegionSpecified ) { return true; }\r
 400 \r
 401         // load next alignment until region overlap is found\r
 402         while ( !IsOverlap(bAlignment) ) {\r
 403             // if no valid alignment available (likely EOF) return failure\r
 404             if ( !LoadNextAlignment(bAlignment) ) { return false; }\r
 405         }\r
 406 \r
 407         // return success (alignment found that overlaps region)\r
 408         return true;\r
 409     }\r
 410 \r
 411     // no valid alignment\r
 412     else { return false; }\r
 413 }\r
 414 \r
 415 // calculate closest indexed file offset for region specified\r
 416 int64_t BamReader::BamReaderPrivate::GetOffset(int refID, int left) {\r
 417 \r
 418     // calculate which bins overlap this region\r
 419     uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);\r
 420     int numBins = BinsFromRegion(refID, left, bins);\r
 421 \r
 422     // get bins for this reference\r
 423     const ReferenceIndex& refIndex = Index.at(refID);\r
 424     const BamBinMap& binMap        = refIndex.Bins;\r
 425 \r
 426     // get minimum offset to consider\r
 427     const LinearOffsetVector& offsets = refIndex.Offsets;\r
 428     uint64_t minOffset = ( (unsigned int)(left>>BAM_LIDX_SHIFT) >= offsets.size() ) ? 0 : offsets.at(left>>BAM_LIDX_SHIFT);\r
 429 \r
 430     // store offsets to beginning of alignment 'chunks'\r
 431     std::vector<int64_t> chunkStarts;\r
 432 \r
 433     // store all alignment 'chunk' starts for bins in this region\r
 434     for (int i = 0; i < numBins; ++i ) {\r
 435         uint16_t binKey = bins[i];\r
 436 \r
 437         map<uint32_t, ChunkVector>::const_iterator binIter = binMap.find(binKey);\r
 438         if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) {\r
 439 \r
 440             const ChunkVector& chunks = (*binIter).second;\r
 441             std::vector<Chunk>::const_iterator chunksIter = chunks.begin();\r
 442             std::vector<Chunk>::const_iterator chunksEnd  = chunks.end();\r
 443             for ( ; chunksIter != chunksEnd; ++chunksIter) {\r
 444                 const Chunk& chunk = (*chunksIter);\r
 445                 if ( chunk.Stop > minOffset ) {\r
 446                     chunkStarts.push_back( chunk.Start );\r
 447                 }\r
 448             }\r
 449         }\r
 450     }\r
 451 \r
 452     // clean up memory\r
 453     free(bins);\r
 454 \r
 455     // if no alignments found, else return smallest offset for alignment starts\r
 456     if ( chunkStarts.size() == 0 ) { return -1; }\r
 457     else { return *min_element(chunkStarts.begin(), chunkStarts.end()); }\r
 458 }\r
 459 \r
 460 // saves BAM bin entry for index\r
 461 void BamReader::BamReaderPrivate::InsertBinEntry(BamBinMap&      binMap,\r
 462                                                  const uint32_t& saveBin,\r
 463                                                  const uint64_t& saveOffset,\r
 464                                                  const uint64_t& lastOffset)\r
 465 {\r
 466     // look up saveBin\r
 467     BamBinMap::iterator binIter = binMap.find(saveBin);\r
 468 \r
 469     // create new chunk\r
 470     Chunk newChunk(saveOffset, lastOffset);\r
 471 \r
 472     // if entry doesn't exist\r
 473     if ( binIter == binMap.end() ) {\r
 474         ChunkVector newChunks;\r
 475         newChunks.push_back(newChunk);\r
 476         binMap.insert( pair<uint32_t, ChunkVector>(saveBin, newChunks));\r
 477     }\r
 478 \r
 479     // otherwise\r
 480     else {\r
 481         ChunkVector& binChunks = (*binIter).second;\r
 482         binChunks.push_back( newChunk );\r
 483     }\r
 484 }\r
 485 \r
 486 // saves linear offset entry for index\r
 487 void BamReader::BamReaderPrivate::InsertLinearOffset(LinearOffsetVector& offsets,\r
 488                                                      const BamAlignment& bAlignment,\r
 489                                                      const uint64_t&     lastOffset)\r
 490 {\r
 491     // get converted offsets\r
 492     int beginOffset = bAlignment.Position >> BAM_LIDX_SHIFT;\r
 493     int endOffset   = ( CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData) - 1) >> BAM_LIDX_SHIFT;\r
 494 \r
 495     // resize vector if necessary\r
 496     int oldSize = offsets.size();\r
 497     int newSize = endOffset + 1;\r
 498     if ( oldSize < newSize ) {        \r
 499         Roundup32(newSize);\r
 500         offsets.resize(newSize, 0);\r
 501     }\r
 502 \r
 503     // store offset\r
 504     for(int i = beginOffset + 1; i <= endOffset ; ++i) {\r
 505         if ( offsets[i] == 0) {\r
 506             offsets[i] = lastOffset;\r
 507         }\r
 508     }\r
 509 }\r
 510 \r
 511 // returns whether alignment overlaps currently specified region (refID, leftBound)\r
 512 bool BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {\r
 513 \r
 514     // if on different reference sequence, quit\r
 515     if ( bAlignment.RefID != CurrentRefID ) { return false; }\r
 516 \r
 517     // read starts after left boundary\r
 518     if ( bAlignment.Position >= CurrentLeft) { return true; }\r
 519 \r
 520     // return whether alignment end overlaps left boundary\r
 521     return ( CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData) >= CurrentLeft );\r
 522 }\r
 523 \r
 524 // jumps to specified region(refID, leftBound) in BAM file, returns success/fail\r
 525 bool BamReader::BamReaderPrivate::Jump(int refID, int position) {\r
 526 \r
 527     // if data exists for this reference and position is valid    \r
 528     if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {\r
 529 \r
 530         // set current region\r
 531         CurrentRefID = refID;\r
 532         CurrentLeft  = position;\r
 533         IsRegionSpecified = true;\r
 534 \r
 535         // calculate offset\r
 536         int64_t offset = GetOffset(CurrentRefID, CurrentLeft);\r
 537 \r
 538         // if in valid offset, return failure\r
 539         if ( offset == -1 ) { return false; }\r
 540 \r
 541         // otherwise return success of seek operation\r
 542         else { return mBGZF.Seek(offset); }\r
 543     }\r
 544 \r
 545     // invalid jump request parameters, return failure\r
 546     return false;\r
 547 }\r
 548 \r
 549 // load BAM header data\r
 550 void BamReader::BamReaderPrivate::LoadHeaderData(void) {\r
 551 \r
 552     // check to see if proper BAM header\r
 553     char buffer[4];\r
 554     if (mBGZF.Read(buffer, 4) != 4) {\r
 555         printf("Could not read header type\n");\r
 556         exit(1);\r
 557     }\r
 558 \r
 559     if (strncmp(buffer, "BAM\001", 4)) {\r
 560         printf("wrong header type!\n");\r
 561         exit(1);\r
 562     }\r
 563 \r
 564     // get BAM header text length\r
 565     mBGZF.Read(buffer, 4);\r
 566     unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);\r
 567     if ( IsBigEndian ) { SwapEndian_32(headerTextLength); }\r
 568     \r
 569     // get BAM header text\r
 570     char* headerText = (char*)calloc(headerTextLength + 1, 1);\r
 571     mBGZF.Read(headerText, headerTextLength);\r
 572     HeaderText = (string)((const char*)headerText);\r
 573 \r
 574     // clean up calloc-ed temp variable\r
 575     free(headerText);\r
 576 }\r
 577 \r
 578 // load existing index data from BAM index file (".bai"), return success/fail\r
 579 bool BamReader::BamReaderPrivate::LoadIndex(void) {\r
 580 \r
 581     // clear out index data\r
 582     ClearIndex();\r
 583 \r
 584     // skip if index file empty\r
 585     if ( IndexFilename.empty() ) { return false; }\r
 586 \r
 587     // open index file, abort on error\r
 588     FILE* indexStream = fopen(IndexFilename.c_str(), "rb");\r
 589     if(!indexStream) {\r
 590         printf("ERROR: Unable to open the BAM index file %s for reading.\n", IndexFilename.c_str() );\r
 591         return false;\r
 592     }\r
 593 \r
 594     size_t elementsRead = 0;\r
 595         \r
 596     // see if index is valid BAM index\r
 597     char magic[4];\r
 598     elementsRead = fread(magic, 1, 4, indexStream);\r
 599     if (strncmp(magic, "BAI\1", 4)) {\r
 600         printf("Problem with index file - invalid format.\n");\r
 601         fclose(indexStream);\r
 602         return false;\r
 603     }\r
 604 \r
 605     // get number of reference sequences\r
 606     uint32_t numRefSeqs;\r
 607     elementsRead = fread(&numRefSeqs, 4, 1, indexStream);\r
 608     if ( IsBigEndian ) { SwapEndian_32(numRefSeqs); }\r
 609     \r
 610     // intialize space for BamIndex data structure\r
 611     Index.reserve(numRefSeqs);\r
 612 \r
 613     // iterate over reference sequences\r
 614     for (unsigned int i = 0; i < numRefSeqs; ++i) {\r
 615 \r
 616         // get number of bins for this reference sequence\r
 617         int32_t numBins;\r
 618         elementsRead = fread(&numBins, 4, 1, indexStream);\r
 619         if ( IsBigEndian ) { SwapEndian_32(numBins); }\r
 620 \r
 621         if (numBins > 0) {\r
 622             RefData& refEntry = References[i];\r
 623             refEntry.RefHasAlignments = true;\r
 624         }\r
 625 \r
 626         // intialize BinVector\r
 627         BamBinMap binMap;\r
 628 \r
 629         // iterate over bins for that reference sequence\r
 630         for (int j = 0; j < numBins; ++j) {\r
 631 \r
 632             // get binID\r
 633             uint32_t binID;\r
 634             elementsRead = fread(&binID, 4, 1, indexStream);\r
 635 \r
 636             // get number of regionChunks in this bin\r
 637             uint32_t numChunks;\r
 638             elementsRead = fread(&numChunks, 4, 1, indexStream);\r
 639 \r
 640             if ( IsBigEndian ) { \r
 641               SwapEndian_32(binID);\r
 642               SwapEndian_32(numChunks);\r
 643             }\r
 644             \r
 645             // intialize ChunkVector\r
 646             ChunkVector regionChunks;\r
 647             regionChunks.reserve(numChunks);\r
 648 \r
 649             // iterate over regionChunks in this bin\r
 650             for (unsigned int k = 0; k < numChunks; ++k) {\r
 651 \r
 652                 // get chunk boundaries (left, right)\r
 653                 uint64_t left;\r
 654                 uint64_t right;\r
 655                 elementsRead = fread(&left, 8, 1, indexStream);\r
 656                 elementsRead = fread(&right, 8, 1, indexStream);\r
 657 \r
 658                 if ( IsBigEndian ) {\r
 659                     SwapEndian_64(left);\r
 660                     SwapEndian_64(right);\r
 661                 }\r
 662                 \r
 663                 // save ChunkPair\r
 664                 regionChunks.push_back( Chunk(left, right) );\r
 665             }\r
 666 \r
 667             // sort chunks for this bin\r
 668             sort( regionChunks.begin(), regionChunks.end(), ChunkLessThan );\r
 669 \r
 670             // save binID, chunkVector for this bin\r
 671             binMap.insert( pair<uint32_t, ChunkVector>(binID, regionChunks) );\r
 672         }\r
 673 \r
 674         // load linear index for this reference sequence\r
 675 \r
 676         // get number of linear offsets\r
 677         int32_t numLinearOffsets;\r
 678         elementsRead = fread(&numLinearOffsets, 4, 1, indexStream);\r
 679         if ( IsBigEndian ) { SwapEndian_32(numLinearOffsets); }\r
 680 \r
 681         // intialize LinearOffsetVector\r
 682         LinearOffsetVector offsets;\r
 683         offsets.reserve(numLinearOffsets);\r
 684 \r
 685         // iterate over linear offsets for this reference sequeence\r
 686         uint64_t linearOffset;\r
 687         for (int j = 0; j < numLinearOffsets; ++j) {\r
 688             // read a linear offset & store\r
 689             elementsRead = fread(&linearOffset, 8, 1, indexStream);\r
 690             if ( IsBigEndian ) { SwapEndian_64(linearOffset); }\r
 691             offsets.push_back(linearOffset);\r
 692         }\r
 693 \r
 694         // sort linear offsets\r
 695         sort( offsets.begin(), offsets.end() );\r
 696 \r
 697         // store index data for that reference sequence\r
 698         Index.push_back( ReferenceIndex(binMap, offsets) );\r
 699     }\r
 700 \r
 701     // close index file (.bai) and return\r
 702     fclose(indexStream);\r
 703     return true;\r
 704 }\r
 705 \r
 706 // populates BamAlignment with alignment data under file pointer, returns success/fail\r
 707 bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {\r
 708 \r
 709     // read in the 'block length' value, make sure it's not zero\r
 710     char buffer[4];\r
 711     mBGZF.Read(buffer, 4);\r
 712     unsigned int blockLength = BgzfData::UnpackUnsignedInt(buffer);\r
 713     if ( IsBigEndian ) { SwapEndian_32(blockLength); }\r
 714     if ( blockLength == 0 ) { return false; }\r
 715 \r
 716     // keep track of bytes read as method progresses\r
 717     int bytesRead = 4;\r
 718 \r
 719     // read in core alignment data, make sure the right size of data was read\r
 720     uint32_t x[8];\r
 721     if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }\r
 722     bytesRead += BAM_CORE_SIZE;\r
 723 \r
 724     if ( IsBigEndian ) {\r
 725         for ( int i = 0; i < 8; ++i ) { \r
 726           SwapEndian_32(x[i]); \r
 727         }\r
 728     }\r
 729     \r
 730     // set BamAlignment 'core' data and character data lengths\r
 731     unsigned int tempValue;\r
 732     unsigned int queryNameLength;\r
 733     unsigned int numCigarOperations;\r
 734     unsigned int querySequenceLength;\r
 735 \r
 736     bAlignment.RefID    = BgzfData::UnpackSignedInt(&x[0]);  \r
 737     bAlignment.Position = BgzfData::UnpackSignedInt(&x[1]);\r
 738     \r
 739     tempValue = BgzfData::UnpackUnsignedInt(&x[2]);\r
 740     bAlignment.Bin        = tempValue >> 16;\r
 741     bAlignment.MapQuality = tempValue >> 8 & 0xff;\r
 742     queryNameLength       = tempValue & 0xff;\r
 743 \r
 744     tempValue = BgzfData::UnpackUnsignedInt(&x[3]);\r
 745     bAlignment.AlignmentFlag = tempValue >> 16;\r
 746     numCigarOperations       = tempValue & 0xffff;\r
 747 \r
 748     querySequenceLength     = BgzfData::UnpackUnsignedInt(&x[4]);\r
 749     bAlignment.MateRefID    = BgzfData::UnpackSignedInt(&x[5]);\r
 750     bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[6]);\r
 751     bAlignment.InsertSize   = BgzfData::UnpackSignedInt(&x[7]);\r
 752     \r
 753     // calculate lengths/offsets\r
 754     const unsigned int dataLength      = blockLength - BAM_CORE_SIZE;\r
 755     const unsigned int cigarDataOffset = queryNameLength;\r
 756     const unsigned int seqDataOffset   = cigarDataOffset + (numCigarOperations * 4);\r
 757     const unsigned int qualDataOffset  = seqDataOffset + (querySequenceLength+1)/2;\r
 758     const unsigned int tagDataOffset   = qualDataOffset + querySequenceLength;\r
 759     const unsigned int tagDataLen      = dataLength - tagDataOffset;\r
 760 \r
 761     // set up destination buffers for character data\r
 762     char* allCharData   = (char*)calloc(sizeof(char), dataLength);\r
 763     uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);\r
 764     char* seqData       = ((char*)allCharData) + seqDataOffset;\r
 765     char* qualData      = ((char*)allCharData) + qualDataOffset;\r
 766     char* tagData       = ((char*)allCharData) + tagDataOffset;\r
 767 \r
 768     // get character data - make sure proper data size was read\r
 769     if ( mBGZF.Read(allCharData, dataLength) != (signed int)dataLength) { return false; }\r
 770     else {\r
 771 \r
 772         bytesRead += dataLength;\r
 773 \r
 774         // clear out any previous string data\r
 775         bAlignment.Name.clear(;)\r
 776         bAlignment.QueryBases.clear();\r
 777         bAlignment.Qualities.clear();\r
 778         bAlignment.AlignedBases.clear();\r
 779         bAlignment.CigarData.clear();\r
 780         bAlignment.TagData.clear();\r
 781 \r
 782         // save name\r
 783         bAlignment.Name = (string)((const char*)(allCharData));\r
 784 \r
 785         // save query sequence\r
 786         // -----------------------\r
 787         // Added: 3-25-2010 DWB\r
 788         // Improved: reduced repeated memory allocations as string grows\r
 789         bAlignment.QueryBases.reserve(querySequenceLength);\r
 790         // -----------------------\r
 791         \r
 792         for (unsigned int i = 0; i < querySequenceLength; ++i) {\r
 793             char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];\r
 794             bAlignment.QueryBases.append( 1, singleBase );\r
 795         }\r
 796 \r
 797         // save sequence length\r
 798         bAlignment.Length = bAlignment.QueryBases.length();\r
 799 \r
 800         // save qualities, convert from numeric QV to FASTQ character\r
 801         // -----------------------\r
 802         // Added: 3-25-2010 DWB\r
 803         // Improved: reduced repeated memory allocations as string grows\r
 804         bAlignment.Qualities.reserve(querySequenceLength);\r
 805         // -----------------------\r
 806         \r
 807         for (unsigned int i = 0; i < querySequenceLength; ++i) {\r
 808             char singleQuality = (char)(qualData[i]+33);\r
 809             bAlignment.Qualities.append( 1, singleQuality );\r
 810         }\r
 811 \r
 812         // save CIGAR-related data;\r
 813         // -----------------------\r
 814         // Added: 3-25-2010 DWB\r
 815         // Improved: reduced repeated memory allocations as string grows\r
 816         bAlignment.AlignedBases.reserve(querySequenceLength);\r
 817         // -----------------------\r
 818         \r
 819         int k = 0;\r
 820         for (unsigned int i = 0; i < numCigarOperations; ++i) {\r
 821 \r
 822             if ( IsBigEndian ) { SwapEndian_32(cigarData[i]); }\r
 823           \r
 824             // build CigarOp struct\r
 825             CigarOp op;\r
 826             op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);\r
 827             op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];\r
 828 \r
 829             // save CigarOp\r
 830             bAlignment.CigarData.push_back(op);\r
 831 \r
 832             // build AlignedBases string\r
 833             switch (op.Type) {\r
 834 \r
 835                 case ('M') :\r
 836                 case ('I') : \r
 837                     bAlignment.AlignedBases.append( bAlignment.QueryBases.substr(k, op.Length) ); // for 'M', 'I' - write bases\r
 838                     // fall through\r
 839                     \r
 840                 case ('S') : \r
 841                     k += op.Length;  // for 'S' - skip over query bases\r
 842                     break;\r
 843 \r
 844                 case ('D') : \r
 845                     bAlignment.AlignedBases.append( op.Length, '-' );   // for 'D' - write gap character\r
 846                     break;\r
 847 \r
 848                 case ('P') : \r
 849                     bAlignment.AlignedBases.append( op.Length, '*' );   // for 'P' - write padding character;\r
 850                     break;\r
 851 \r
 852                 case ('N') : \r
 853                     bAlignment.AlignedBases.append( op.Length, 'N' );  // for 'N' - write N's, skip bases in query sequence\r
 854                     // -----------------------\r
 855                     // Removed: 3-25-2010 DWB\r
 856                     // Contributed: ARQ\r
 857                     // Fixed: compliance with actual 'N' definition in BAM spec\r
 858                     // k += op.Length;\r
 859                     // -----------------------\r
 860                     break;\r
 861 \r
 862                 case ('H') : \r
 863                     break;  // for 'H' - do nothing, move to next op\r
 864 \r
 865                 default : \r
 866                     printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here\r
 867                     free(allCharData);\r
 868                     exit(1);\r
 869             }\r
 870         }\r
 871 \r
 872         // -----------------------\r
 873         // Added: 3-25-2010 DWB\r
 874         // Fixed: endian-correctness for tag data\r
 875         // -----------------------\r
 876         if ( IsBigEndian ) {\r
 877             int i = 0;\r
 878             while ( i < tagDataLen ) {\r
 879                 \r
 880                 i += 2;                                 // skip tag type (e.g. "RG", "NM", etc)\r
 881                 uint8_t type = toupper(tagData[i]);     // lower & upper case letters have same meaning \r
 882                 ++i;                                    // skip value type\r
 883                 \r
 884                 switch (type) {\r
 885                   \r
 886                     case('A') :\r
 887                     case('C') : \r
 888                         ++i;\r
 889                         break;\r
 890                                 \r
 891                     case('S') : \r
 892                         SwapEndian_16p(&tagData[i]); \r
 893                         i+=2; // sizeof(uint16_t)\r
 894                         break;\r
 895                                 \r
 896                     case('F') :\r
 897                     case('I') : \r
 898                         SwapEndian_32p(&tagData[i]);\r
 899                         i+=4; // sizeof(uint32_t)\r
 900                         break;\r
 901                                 \r
 902                     case('D') : \r
 903                         SwapEndian_64p(&tagData[i]);\r
 904                         i+=8; // sizeof(uint64_t)\r
 905                         break;\r
 906                                 \r
 907                     case('H') :\r
 908                     case('Z') : \r
 909                         while (tagData[i]) { ++i; }\r
 910                         ++i; // increment one more for null terminator\r
 911                         break;\r
 912                                 \r
 913                     default : \r
 914                         printf("ERROR: Invalid tag value type\n"); // shouldn't get here\r
 915                         free(allCharData);\r
 916                         exit(1); \r
 917                 }\r
 918             }\r
 919         }\r
 920         \r
 921         // store tag data\r
 922         bAlignment.TagData.resize(tagDataLen);\r
 923         memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLen);\r
 924     }\r
 925 \r
 926     free(allCharData);\r
 927     return true;\r
 928 }\r
 929 \r
 930 // loads reference data from BAM file\r
 931 void BamReader::BamReaderPrivate::LoadReferenceData(void) {\r
 932 \r
 933     // get number of reference sequences\r
 934     char buffer[4];\r
 935     mBGZF.Read(buffer, 4);\r
 936     unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);\r
 937     if ( IsBigEndian ) { SwapEndian_32(numberRefSeqs); }\r
 938     if (numberRefSeqs == 0) { return; }\r
 939     References.reserve((int)numberRefSeqs);\r
 940 \r
 941     // iterate over all references in header\r
 942     for (unsigned int i = 0; i != numberRefSeqs; ++i) {\r
 943 \r
 944         // get length of reference name\r
 945         mBGZF.Read(buffer, 4);\r
 946         unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);\r
 947         if ( IsBigEndian ) { SwapEndian_32(refNameLength); }\r
 948         char* refName = (char*)calloc(refNameLength, 1);\r
 949 \r
 950         // get reference name and reference sequence length\r
 951         mBGZF.Read(refName, refNameLength);\r
 952         mBGZF.Read(buffer, 4);\r
 953         int refLength = BgzfData::UnpackSignedInt(buffer);\r
 954         if ( IsBigEndian ) { SwapEndian_32(refLength); }\r
 955 \r
 956         // store data for reference\r
 957         RefData aReference;\r
 958         aReference.RefName   = (string)((const char*)refName);\r
 959         aReference.RefLength = refLength;\r
 960         References.push_back(aReference);\r
 961 \r
 962         // clean up calloc-ed temp variable\r
 963         free(refName);\r
 964     }\r
 965 }\r
 966 \r
 967 // merges 'alignment chunks' in BAM bin (used for index building)\r
 968 void BamReader::BamReaderPrivate::MergeChunks(void) {\r
 969 \r
 970     // iterate over reference enties\r
 971     BamIndex::iterator indexIter = Index.begin();\r
 972     BamIndex::iterator indexEnd  = Index.end();\r
 973     for ( ; indexIter != indexEnd; ++indexIter ) {\r
 974 \r
 975         // get BAM bin map for this reference\r
 976         ReferenceIndex& refIndex = (*indexIter);\r
 977         BamBinMap& bamBinMap = refIndex.Bins;\r
 978 \r
 979         // iterate over BAM bins\r
 980         BamBinMap::iterator binIter = bamBinMap.begin();\r
 981         BamBinMap::iterator binEnd  = bamBinMap.end();\r
 982         for ( ; binIter != binEnd; ++binIter ) {\r
 983 \r
 984             // get chunk vector for this bin\r
 985             ChunkVector& binChunks = (*binIter).second;\r
 986             if ( binChunks.size() == 0 ) { continue; }\r
 987 \r
 988             ChunkVector mergedChunks;\r
 989             mergedChunks.push_back( binChunks[0] );\r
 990 \r
 991             // iterate over chunks\r
 992             int i = 0;\r
 993             ChunkVector::iterator chunkIter = binChunks.begin();\r
 994             ChunkVector::iterator chunkEnd  = binChunks.end();\r
 995             for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {\r
 996 \r
 997                 // get 'currentChunk' based on numeric index\r
 998                 Chunk& currentChunk = mergedChunks[i];\r
 999 \r
1000                 // get iteratorChunk based on vector iterator\r
1001                 Chunk& iteratorChunk = (*chunkIter);\r
1002 \r
1003                 // if currentChunk.Stop(shifted) == iterator Chunk.Start(shifted)\r
1004                 if ( currentChunk.Stop>>16 == iteratorChunk.Start>>16 ) {\r
1005 \r
1006                     // set currentChunk.Stop to iteratorChunk.Stop\r
1007                     currentChunk.Stop = iteratorChunk.Stop;\r
1008                 }\r
1009 \r
1010                 // otherwise\r
1011                 else {\r
1012                     // set currentChunk + 1 to iteratorChunk\r
1013                     mergedChunks.push_back(iteratorChunk);\r
1014                     ++i;\r
1015                 }\r
1016             }\r
1017 \r
1018             // saved merged chunk vector\r
1019             (*binIter).second = mergedChunks;\r
1020         }\r
1021     }\r
1022 }\r
1023 \r
1024 // opens BAM file (and index)\r
1025 void BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename) {\r
1026 \r
1027     Filename = filename;\r
1028     IndexFilename = indexFilename;\r
1029 \r
1030     // open the BGZF file for reading, retrieve header text & reference data\r
1031     mBGZF.Open(filename, "rb");\r
1032     LoadHeaderData();\r
1033     LoadReferenceData();\r
1034 \r
1035     // store file offset of first alignment\r
1036     AlignmentsBeginOffset = mBGZF.Tell();\r
1037 \r
1038     // open index file & load index data (if exists)\r
1039     if ( !IndexFilename.empty() ) {\r
1040         LoadIndex();\r
1041     }\r
1042 }\r
1043 \r
1044 // returns BAM file pointer to beginning of alignment data\r
1045 bool BamReader::BamReaderPrivate::Rewind(void) {\r
1046 \r
1047     // find first reference that has alignments in the BAM file\r
1048     int refID = 0;\r
1049     int refCount = References.size();\r
1050     for ( ; refID < refCount; ++refID ) {\r
1051         if ( References.at(refID).RefHasAlignments ) { break; }\r
1052     }\r
1053 \r
1054     // store default bounds for first alignment\r
1055     CurrentRefID = refID;\r
1056     CurrentLeft = 0;\r
1057     IsRegionSpecified = false;\r
1058 \r
1059     // return success/failure of seek\r
1060     return mBGZF.Seek(AlignmentsBeginOffset);\r
1061 }\r
1062 \r
1063 // rounds value up to next power-of-2 (used in index building)\r
1064 void BamReader::BamReaderPrivate::Roundup32(int& value) {    \r
1065     --value;\r
1066     value |= value >> 1;\r
1067     value |= value >> 2;\r
1068     value |= value >> 4;\r
1069     value |= value >> 8;\r
1070     value |= value >> 16;\r
1071     ++value;\r
1072 }\r
1073 \r
1074 // saves index data to BAM index file (".bai"), returns success/fail\r
1075 bool BamReader::BamReaderPrivate::WriteIndex(void) {\r
1076 \r
1077     IndexFilename = Filename + ".bai";\r
1078     FILE* indexStream = fopen(IndexFilename.c_str(), "wb");\r
1079     if ( indexStream == 0 ) {\r
1080         printf("ERROR: Could not open file to save index\n");\r
1081         return false;\r
1082     }\r
1083 \r
1084     // write BAM index header\r
1085     fwrite("BAI\1", 1, 4, indexStream);\r
1086 \r
1087     // write number of reference sequences\r
1088     int32_t numReferenceSeqs = Index.size();\r
1089     if ( IsBigEndian ) { SwapEndian_32(numReferenceSeqs); }\r
1090     fwrite(&numReferenceSeqs, 4, 1, indexStream);\r
1091 \r
1092     // iterate over reference sequences\r
1093     BamIndex::const_iterator indexIter = Index.begin();\r
1094     BamIndex::const_iterator indexEnd  = Index.end();\r
1095     for ( ; indexIter != indexEnd; ++ indexIter ) {\r
1096 \r
1097         // get reference index data\r
1098         const ReferenceIndex& refIndex = (*indexIter);\r
1099         const BamBinMap& binMap = refIndex.Bins;\r
1100         const LinearOffsetVector& offsets = refIndex.Offsets;\r
1101 \r
1102         // write number of bins\r
1103         int32_t binCount = binMap.size();\r
1104         if ( IsBigEndian ) { SwapEndian_32(binCount); }\r
1105         fwrite(&binCount, 4, 1, indexStream);\r
1106 \r
1107         // iterate over bins\r
1108         BamBinMap::const_iterator binIter = binMap.begin();\r
1109         BamBinMap::const_iterator binEnd  = binMap.end();\r
1110         for ( ; binIter != binEnd; ++binIter ) {\r
1111 \r
1112             // get bin data (key and chunk vector)\r
1113             uint32_t binKey = (*binIter).first;\r
1114             const ChunkVector& binChunks = (*binIter).second;\r
1115 \r
1116             // save BAM bin key\r
1117             if ( IsBigEndian ) { SwapEndian_32(binKey); }\r
1118             fwrite(&binKey, 4, 1, indexStream);\r
1119 \r
1120             // save chunk count\r
1121             int32_t chunkCount = binChunks.size();\r
1122             if ( IsBigEndian ) { SwapEndian_32(chunkCount); }\r
1123             fwrite(&chunkCount, 4, 1, indexStream);\r
1124 \r
1125             // iterate over chunks\r
1126             ChunkVector::const_iterator chunkIter = binChunks.begin();\r
1127             ChunkVector::const_iterator chunkEnd  = binChunks.end();\r
1128             for ( ; chunkIter != chunkEnd; ++chunkIter ) {\r
1129 \r
1130                 // get current chunk data\r
1131                 const Chunk& chunk    = (*chunkIter);\r
1132                 uint64_t start = chunk.Start;\r
1133                 uint64_t stop  = chunk.Stop;\r
1134 \r
1135                 if ( IsBigEndian ) {\r
1136                     SwapEndian_64(start);\r
1137                     SwapEndian_64(stop);\r
1138                 }\r
1139                 \r
1140                 // save chunk offsets\r
1141                 fwrite(&start, 8, 1, indexStream);\r
1142                 fwrite(&stop,  8, 1, indexStream);\r
1143             }\r
1144         }\r
1145 \r
1146         // write linear offsets size\r
1147         int32_t offsetSize = offsets.size();\r
1148         if ( IsBigEndian ) { SwapEndian_32(offsetSize); }\r
1149         fwrite(&offsetSize, 4, 1, indexStream);\r
1150 \r
1151         // iterate over linear offsets\r
1152         LinearOffsetVector::const_iterator offsetIter = offsets.begin();\r
1153         LinearOffsetVector::const_iterator offsetEnd  = offsets.end();\r
1154         for ( ; offsetIter != offsetEnd; ++offsetIter ) {\r
1155 \r
1156             // write linear offset value\r
1157             uint64_t linearOffset = (*offsetIter);\r
1158             if ( IsBigEndian ) { SwapEndian_64(linearOffset); }\r
1159             fwrite(&linearOffset, 8, 1, indexStream);\r
1160         }\r
1161     }\r
1162 \r
1163     // flush buffer, close file, and return success\r
1164     fflush(indexStream);\r
1165     fclose(indexStream);\r
1166     return true;\r
1167 }\r