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