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