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