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