src/api/BamReader.cpp

   1 // ***************************************************************************\r
   2 // BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg\r
   3 // Marth Lab, Department of Biology, Boston College\r
   4 // All rights reserved.\r
   5 // ---------------------------------------------------------------------------\r
   6 // Last modified: 30 August 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 #include <iostream>\r
  20 \r
  21 // BamTools includes\r
  22 #include "BGZF.h"\r
  23 #include "BamReader.h"\r
  24 #include "BamIndex.h"\r
  25 using namespace BamTools;\r
  26 using namespace std;\r
  27 \r
  28 struct BamReader::BamReaderPrivate {\r
  29 \r
  30     // -------------------------------\r
  31     // structs, enums, typedefs\r
  32     // -------------------------------\r
  33     enum RegionState { BEFORE_REGION = 0\r
  34                      , WITHIN_REGION\r
  35                      , AFTER_REGION\r
  36                      };\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     BamIndex* NewIndex;\r
  47     RefVector References;\r
  48     bool      IsIndexLoaded;\r
  49     int64_t   AlignmentsBeginOffset;\r
  50     string    Filename;\r
  51     string    IndexFilename;\r
  52     \r
  53     // system data\r
  54     bool IsBigEndian;\r
  55 \r
  56     // user-specified region values\r
  57     BamRegion Region;\r
  58     bool IsLeftBoundSpecified;\r
  59     bool IsRightBoundSpecified;\r
  60     \r
  61     bool IsRegionSpecified;\r
  62     int  CurrentRefID;\r
  63     int  CurrentLeft;\r
  64 \r
  65     // parent BamReader\r
  66     BamReader* Parent;\r
  67     \r
  68     // BAM character constants\r
  69     const char* DNA_LOOKUP;\r
  70     const char* CIGAR_LOOKUP;\r
  71 \r
  72     // -------------------------------\r
  73     // constructor & destructor\r
  74     // -------------------------------\r
  75     BamReaderPrivate(BamReader* parent);\r
  76     ~BamReaderPrivate(void);\r
  77 \r
  78     // -------------------------------\r
  79     // "public" interface\r
  80     // -------------------------------\r
  81 \r
  82     // file operations\r
  83     void Close(void);\r
  84     bool Jump(int refID, int position = 0);\r
  85     bool Open(const string& filename, const string& indexFilename = "");\r
  86     bool Rewind(void);\r
  87     bool SetRegion(const BamRegion& region);\r
  88 \r
  89     // access alignment data\r
  90     bool GetNextAlignment(BamAlignment& bAlignment);\r
  91     bool GetNextAlignmentCore(BamAlignment& bAlignment);\r
  92 \r
  93     // access auxiliary data\r
  94     int GetReferenceID(const string& refName) const;\r
  95 \r
  96     // index operations\r
  97     bool CreateIndex(bool useDefaultIndex);\r
  98 \r
  99     // -------------------------------\r
 100     // internal methods\r
 101     // -------------------------------\r
 102 \r
 103     // *** reading alignments and auxiliary data *** //\r
 104 \r
 105     // fills out character data for BamAlignment data\r
 106     bool BuildCharData(BamAlignment& bAlignment);\r
 107     // checks to see if alignment overlaps current region\r
 108     RegionState IsOverlap(BamAlignment& bAlignment);\r
 109     // retrieves header text from BAM file\r
 110     void LoadHeaderData(void);\r
 111     // retrieves BAM alignment under file pointer\r
 112     bool LoadNextAlignment(BamAlignment& bAlignment);\r
 113     // builds reference data structure from BAM file\r
 114     void LoadReferenceData(void);\r
 115 \r
 116     // *** index file handling *** //\r
 117 \r
 118     // clear out inernal index data structure\r
 119     void ClearIndex(void);\r
 120     // loads index from BAM index file\r
 121     bool LoadIndex(void);\r
 122 };\r
 123 \r
 124 // -----------------------------------------------------\r
 125 // BamReader implementation (wrapper around BRPrivate)\r
 126 // -----------------------------------------------------\r
 127 // constructor\r
 128 BamReader::BamReader(void) {\r
 129     d = new BamReaderPrivate(this);\r
 130 }\r
 131 \r
 132 // destructor\r
 133 BamReader::~BamReader(void) {\r
 134     delete d;\r
 135     d = 0;\r
 136 }\r
 137 \r
 138 // file operations\r
 139 void BamReader::Close(void) { d->Close(); }\r
 140 bool BamReader::IsIndexLoaded(void) const { return d->IsIndexLoaded; }\r
 141 bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }\r
 142 bool BamReader::Jump(int refID, int position) { \r
 143     d->Region.LeftRefID = refID;\r
 144     d->Region.LeftPosition = position;\r
 145     d->IsLeftBoundSpecified = true;\r
 146     d->IsRightBoundSpecified = false;\r
 147     return d->Jump(refID, position); \r
 148 }\r
 149 bool BamReader::Open(const string& filename, const string& indexFilename) { return d->Open(filename, indexFilename); }\r
 150 bool BamReader::Rewind(void) { return d->Rewind(); }\r
 151 bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }\r
 152 bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {\r
 153     return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );\r
 154 }\r
 155 \r
 156 // access alignment data\r
 157 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }\r
 158 bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }\r
 159 \r
 160 // access auxiliary data\r
 161 const string BamReader::GetHeaderText(void) const { return d->HeaderText; }\r
 162 int BamReader::GetReferenceCount(void) const { return d->References.size(); }\r
 163 const RefVector& BamReader::GetReferenceData(void) const { return d->References; }\r
 164 int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }\r
 165 const std::string BamReader::GetFilename(void) const { return d->Filename; }\r
 166 \r
 167 // index operations\r
 168 bool BamReader::CreateIndex(bool useDefaultIndex) { return d->CreateIndex(useDefaultIndex); }\r
 169 \r
 170 // -----------------------------------------------------\r
 171 // BamReaderPrivate implementation\r
 172 // -----------------------------------------------------\r
 173 \r
 174 // constructor\r
 175 BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)\r
 176     : NewIndex(0)\r
 177     , IsIndexLoaded(false)\r
 178     , AlignmentsBeginOffset(0)\r
 179     , IsLeftBoundSpecified(false)\r
 180     , IsRightBoundSpecified(false)\r
 181     , IsRegionSpecified(false)\r
 182     , CurrentRefID(0)\r
 183     , CurrentLeft(0)\r
 184     , Parent(parent)\r
 185     , DNA_LOOKUP("=ACMGRSVTWYHKDBN")\r
 186     , CIGAR_LOOKUP("MIDNSHP")\r
 187 { \r
 188     IsBigEndian = SystemIsBigEndian();\r
 189 }\r
 190 \r
 191 // destructor\r
 192 BamReader::BamReaderPrivate::~BamReaderPrivate(void) {\r
 193     Close();\r
 194 }\r
 195 \r
 196 bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {\r
 197   \r
 198     // calculate character lengths/offsets\r
 199     const unsigned int dataLength      = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;\r
 200 //     const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;\r
 201     const unsigned int seqDataOffset   = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);\r
 202     const unsigned int qualDataOffset  = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;\r
 203     const unsigned int tagDataOffset   = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;\r
 204     const unsigned int tagDataLength   = dataLength - tagDataOffset;\r
 205       \r
 206     // set up char buffers\r
 207     const char*     allCharData = bAlignment.SupportData.AllCharData.data();\r
 208 //           uint32_t* cigarData   = (uint32_t*)(allCharData + cigarDataOffset);\r
 209     const char*     seqData     = ((const char*)allCharData) + seqDataOffset;\r
 210     const char*     qualData    = ((const char*)allCharData) + qualDataOffset;\r
 211           char*     tagData     = ((char*)allCharData) + tagDataOffset;\r
 212   \r
 213     // store alignment name (relies on null char in name as terminator)\r
 214     bAlignment.Name.assign((const char*)(allCharData));    \r
 215           \r
 216 //     // save CigarOps \r
 217 //     CigarOp op;\r
 218 //     bAlignment.CigarData.clear();\r
 219 //     bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);\r
 220 //     for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {\r
 221 // \r
 222 //         // swap if necessary\r
 223 //         if ( IsBigEndian ) { SwapEndian_32(cigarData[i]); }\r
 224 //       \r
 225 //         // build CigarOp structure\r
 226 //         op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);\r
 227 //         op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];\r
 228 // \r
 229 //         // save CigarOp\r
 230 //         bAlignment.CigarData.push_back(op);\r
 231 //     }\r
 232           \r
 233     // save query sequence\r
 234     bAlignment.QueryBases.clear();\r
 235     bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);\r
 236     for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {\r
 237         char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];\r
 238         bAlignment.QueryBases.append(1, singleBase);\r
 239     }\r
 240   \r
 241     // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character\r
 242     bAlignment.Qualities.clear();\r
 243     bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);\r
 244     for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {\r
 245         char singleQuality = (char)(qualData[i]+33);\r
 246         bAlignment.Qualities.append(1, singleQuality);\r
 247     }\r
 248     \r
 249     // if QueryBases is empty (and this is a allowed case)\r
 250     if ( bAlignment.QueryBases.empty() ) \r
 251         bAlignment.AlignedBases = bAlignment.QueryBases;\r
 252     \r
 253     // if QueryBases contains data, then build AlignedBases using CIGAR data\r
 254     else {\r
 255     \r
 256         // resize AlignedBases\r
 257         bAlignment.AlignedBases.clear();\r
 258         bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);\r
 259       \r
 260         // iterate over CigarOps\r
 261         int k = 0;\r
 262         vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();\r
 263         vector<CigarOp>::const_iterator cigarEnd  = bAlignment.CigarData.end();\r
 264         for ( ; cigarIter != cigarEnd; ++cigarIter ) {\r
 265             \r
 266             const CigarOp& op = (*cigarIter);\r
 267             switch(op.Type) {\r
 268               \r
 269                 case ('M') :\r
 270                 case ('I') :\r
 271                     bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases\r
 272                     // fall through\r
 273                 \r
 274                 case ('S') :\r
 275                     k += op.Length;                                     // for 'S' - soft clip, skip over query bases\r
 276                     break;\r
 277                     \r
 278                 case ('D') :\r
 279                     bAlignment.AlignedBases.append(op.Length, '-');     // for 'D' - write gap character\r
 280                     break;\r
 281                     \r
 282                 case ('P') :\r
 283                     bAlignment.AlignedBases.append( op.Length, '*' );   // for 'P' - write padding character\r
 284                     break;\r
 285                     \r
 286                 case ('N') :\r
 287                     bAlignment.AlignedBases.append( op.Length, 'N' );  // for 'N' - write N's, skip bases in original query sequence\r
 288                     break;\r
 289                     \r
 290                 case ('H') :\r
 291                     break;  // for 'H' - hard clip, do nothing to AlignedBases, move to next op\r
 292                     \r
 293                 default:\r
 294                     printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here\r
 295                     exit(1);\r
 296             }\r
 297         }\r
 298     }\r
 299  \r
 300     // -----------------------\r
 301     // Added: 3-25-2010 DB\r
 302     // Fixed: endian-correctness for tag data\r
 303     // -----------------------\r
 304     if ( IsBigEndian ) {\r
 305         int i = 0;\r
 306         while ( (unsigned int)i < tagDataLength ) {\r
 307           \r
 308             i += 2; // skip tag type (e.g. "RG", "NM", etc)\r
 309             uint8_t type = toupper(tagData[i]);     // lower & upper case letters have same meaning \r
 310             ++i;                                    // skip value type\r
 311     \r
 312             switch (type) {\r
 313                 \r
 314                 case('A') :\r
 315                 case('C') : \r
 316                     ++i;\r
 317                     break;\r
 318 \r
 319                 case('S') : \r
 320                     SwapEndian_16p(&tagData[i]); \r
 321                     i += sizeof(uint16_t);\r
 322                     break;\r
 323                     \r
 324                 case('F') :\r
 325                 case('I') : \r
 326                     SwapEndian_32p(&tagData[i]);\r
 327                     i += sizeof(uint32_t);\r
 328                     break;\r
 329                 \r
 330                 case('D') : \r
 331                     SwapEndian_64p(&tagData[i]);\r
 332                     i += sizeof(uint64_t);\r
 333                     break;\r
 334                 \r
 335                 case('H') :\r
 336                 case('Z') : \r
 337                     while (tagData[i]) { ++i; }\r
 338                     ++i; // increment one more for null terminator\r
 339                     break;\r
 340                 \r
 341                 default : \r
 342                     printf("ERROR: Invalid tag value type\n"); // shouldn't get here\r
 343                     exit(1);\r
 344             }\r
 345         }\r
 346     }\r
 347     \r
 348     // store TagData\r
 349     bAlignment.TagData.clear();\r
 350     bAlignment.TagData.resize(tagDataLength);\r
 351     memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);\r
 352     \r
 353     // clear the core-only flag\r
 354     bAlignment.SupportData.HasCoreOnly = false;\r
 355     \r
 356     // return success\r
 357     return true;\r
 358 }\r
 359 \r
 360 // clear index data structure\r
 361 void BamReader::BamReaderPrivate::ClearIndex(void) {\r
 362     delete NewIndex;\r
 363     NewIndex = 0;\r
 364     IsIndexLoaded = false;\r
 365 }\r
 366 \r
 367 // closes the BAM file\r
 368 void BamReader::BamReaderPrivate::Close(void) {\r
 369     \r
 370     // close BGZF file stream\r
 371     mBGZF.Close();\r
 372     \r
 373     // clear out index data\r
 374     ClearIndex();\r
 375     \r
 376     // clear out header data\r
 377     HeaderText.clear();\r
 378     \r
 379     // clear out region flags\r
 380     IsLeftBoundSpecified  = false;\r
 381     IsRightBoundSpecified = false;\r
 382     IsRegionSpecified     = false;\r
 383 }\r
 384 \r
 385 // create BAM index from BAM file (keep structure in memory) and write to default index output file\r
 386 bool BamReader::BamReaderPrivate::CreateIndex(bool useDefaultIndex) {\r
 387 \r
 388     // clear out prior index data\r
 389     ClearIndex();\r
 390     \r
 391     // create default index\r
 392     if ( useDefaultIndex )\r
 393         NewIndex = new BamDefaultIndex(&mBGZF, Parent, IsBigEndian);\r
 394     // create BamTools 'custom' index\r
 395     else\r
 396         NewIndex = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);\r
 397     \r
 398     // build new index\r
 399     bool ok = true;\r
 400     ok &= NewIndex->Build();\r
 401     IsIndexLoaded = ok;\r
 402     \r
 403     // attempt to save index data to file\r
 404     ok &= NewIndex->Write(Filename); \r
 405     \r
 406     // return success/fail of both building & writing index\r
 407     return ok;\r
 408 }\r
 409 \r
 410 // get next alignment (from specified region, if given)\r
 411 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {\r
 412 \r
 413     // if valid alignment found, attempt to parse char data, and return success/failure\r
 414     if ( GetNextAlignmentCore(bAlignment) )\r
 415         return BuildCharData(bAlignment);\r
 416     \r
 417     // no valid alignment found\r
 418     else\r
 419         return false;\r
 420 }\r
 421 \r
 422 // retrieves next available alignment core data (returns success/fail)\r
 423 // ** DOES NOT parse any character data (read name, bases, qualities, tag data)\r
 424 //    these can be accessed, if necessary, from the supportData \r
 425 // useful for operations requiring ONLY positional or other alignment-related information\r
 426 bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {\r
 427 \r
 428     // if valid alignment available\r
 429     if ( LoadNextAlignment(bAlignment) ) {\r
 430 \r
 431         // set core-only flag\r
 432         bAlignment.SupportData.HasCoreOnly = true;\r
 433       \r
 434         // if region not specified, return success\r
 435         if ( !IsLeftBoundSpecified ) return true;\r
 436 \r
 437         // determine region state (before, within, after)\r
 438         BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);\r
 439       \r
 440         // if alignment lies after region, return false\r
 441         if ( state == AFTER_REGION ) \r
 442             return false;\r
 443 \r
 444         while ( state != WITHIN_REGION ) {\r
 445             // if no valid alignment available (likely EOF) return failure\r
 446             if ( !LoadNextAlignment(bAlignment) ) return false;\r
 447             // if alignment lies after region, return false (no available read within region)\r
 448             state = IsOverlap(bAlignment);\r
 449             if ( state == AFTER_REGION) return false;\r
 450             \r
 451         }\r
 452 \r
 453         // return success (alignment found that overlaps region)\r
 454         return true;\r
 455     }\r
 456 \r
 457     // no valid alignment\r
 458     else\r
 459         return false;\r
 460 }\r
 461 \r
 462 // returns RefID for given RefName (returns References.size() if not found)\r
 463 int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {\r
 464 \r
 465     // retrieve names from reference data\r
 466     vector<string> refNames;\r
 467     RefVector::const_iterator refIter = References.begin();\r
 468     RefVector::const_iterator refEnd  = References.end();\r
 469     for ( ; refIter != refEnd; ++refIter) {\r
 470         refNames.push_back( (*refIter).RefName );\r
 471     }\r
 472 \r
 473     // return 'index-of' refName ( if not found, returns refNames.size() )\r
 474     return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));\r
 475 }\r
 476 \r
 477 // returns region state - whether alignment ends before, overlaps, or starts after currently specified region\r
 478 // this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true\r
 479 BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {\r
 480     \r
 481     // --------------------------------------------------\r
 482     // check alignment start against right bound cutoff\r
 483     \r
 484     // if full region of interest was given\r
 485     if ( IsRightBoundSpecified ) {\r
 486       \r
 487         // read starts on right bound reference, but AFTER right bound position\r
 488         if ( bAlignment.RefID == Region.RightRefID && bAlignment.Position > Region.RightPosition )\r
 489             return AFTER_REGION;\r
 490       \r
 491         // if read starts on reference AFTER right bound, return false\r
 492         if ( bAlignment.RefID > Region.RightRefID ) \r
 493             return AFTER_REGION;\r
 494     }\r
 495   \r
 496     // --------------------------------------------------------\r
 497     // no right bound given OR read starts before right bound\r
 498     // so, check if it overlaps left bound \r
 499   \r
 500     // if read starts on left bound reference AND after left boundary, return success\r
 501     if ( bAlignment.RefID == Region.LeftRefID && bAlignment.Position >= Region.LeftPosition)\r
 502         return WITHIN_REGION;\r
 503   \r
 504     // if read is on any reference sequence before left bound, return false\r
 505     if ( bAlignment.RefID < Region.LeftRefID )\r
 506         return BEFORE_REGION;\r
 507 \r
 508     // --------------------------------------------------------\r
 509     // read is on left bound reference, but starts before left bound position\r
 510 \r
 511     // if it overlaps, return WITHIN_REGION\r
 512     if ( bAlignment.GetEndPosition() >= Region.LeftPosition )\r
 513         return WITHIN_REGION;\r
 514     // else begins before left bound position\r
 515     else\r
 516         return BEFORE_REGION;\r
 517 }\r
 518 \r
 519 // jumps to specified region(refID, leftBound) in BAM file, returns success/fail\r
 520 bool BamReader::BamReaderPrivate::Jump(int refID, int position) {\r
 521 \r
 522     // -----------------------------------------------------------------------\r
 523     // check for existing index \r
 524     if ( !IsIndexLoaded || NewIndex == 0 ) return false; \r
 525     // see if reference has alignments\r
 526     if ( !NewIndex->HasAlignments(refID) ) return false; \r
 527     // make sure position is valid\r
 528     if ( position > References.at(refID).RefLength ) return false;\r
 529     \r
 530     // determine possible offsets\r
 531     vector<int64_t> offsets;\r
 532     if ( !NewIndex->GetOffsets(Region, IsRightBoundSpecified, offsets) ) {\r
 533         printf("ERROR: Could not jump: unable to calculate offset for specified region.\n");\r
 534         return false;\r
 535     }\r
 536       \r
 537     // iterate through offsets\r
 538     BamAlignment bAlignment;\r
 539     bool result = true;\r
 540     for ( vector<int64_t>::const_iterator o = offsets.begin(); o != offsets.end(); ++o) {\r
 541         \r
 542         // attempt seek & load first available alignment\r
 543         result &= mBGZF.Seek(*o);\r
 544         LoadNextAlignment(bAlignment);\r
 545         \r
 546         // if this alignment corresponds to desired position\r
 547         // return success of seeking back to 'current offset'\r
 548         if ( (bAlignment.RefID == refID && bAlignment.Position + bAlignment.Length > position) || (bAlignment.RefID > refID) ) {\r
 549             if ( o != offsets.begin() ) --o;\r
 550             return mBGZF.Seek(*o);\r
 551         }\r
 552     }\r
 553     \r
 554     return result;\r
 555 }\r
 556 \r
 557 // load BAM header data\r
 558 void BamReader::BamReaderPrivate::LoadHeaderData(void) {\r
 559 \r
 560     // check to see if proper BAM header\r
 561     char buffer[4];\r
 562     if (mBGZF.Read(buffer, 4) != 4) {\r
 563         printf("Could not read header type\n");\r
 564         exit(1);\r
 565     }\r
 566 \r
 567     if (strncmp(buffer, "BAM\001", 4)) {\r
 568         printf("wrong header type!\n");\r
 569         exit(1);\r
 570     }\r
 571 \r
 572     // get BAM header text length\r
 573     mBGZF.Read(buffer, 4);\r
 574     unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);\r
 575     if ( IsBigEndian ) { SwapEndian_32(headerTextLength); }\r
 576     \r
 577     // get BAM header text\r
 578     char* headerText = (char*)calloc(headerTextLength + 1, 1);\r
 579     mBGZF.Read(headerText, headerTextLength);\r
 580     HeaderText = (string)((const char*)headerText);\r
 581 \r
 582     // clean up calloc-ed temp variable\r
 583     free(headerText);\r
 584 }\r
 585 \r
 586 // load existing index data from BAM index file (".bai"), return success/fail\r
 587 bool BamReader::BamReaderPrivate::LoadIndex(void) {\r
 588 \r
 589     // clear out any existing index data\r
 590     ClearIndex();\r
 591 \r
 592     // skip if index file empty\r
 593     if ( IndexFilename.empty() )\r
 594         return false;\r
 595 \r
 596     // check supplied filename for index type\r
 597     size_t defaultExtensionFound = IndexFilename.find(".bai");\r
 598     size_t customExtensionFound  = IndexFilename.find(".bti");\r
 599     \r
 600     // if SAM/BAM default (".bai")\r
 601     if ( defaultExtensionFound != string::npos )\r
 602         NewIndex = new BamDefaultIndex(&mBGZF, Parent, IsBigEndian);\r
 603     \r
 604     // if BamTools custom index (".bti")\r
 605     else if ( customExtensionFound != string::npos )\r
 606         NewIndex = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);\r
 607     \r
 608     // else unknown\r
 609     else {\r
 610         printf("ERROR: Unknown index file extension.\n");\r
 611         return false;\r
 612     }\r
 613     \r
 614     // return success of loading index data from file\r
 615     IsIndexLoaded = NewIndex->Load(IndexFilename);\r
 616     return IsIndexLoaded;\r
 617 }\r
 618 \r
 619 // populates BamAlignment with alignment data under file pointer, returns success/fail\r
 620 bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {\r
 621 \r
 622     // read in the 'block length' value, make sure it's not zero\r
 623     char buffer[4];\r
 624     mBGZF.Read(buffer, 4);\r
 625     bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);\r
 626     if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }\r
 627     if ( bAlignment.SupportData.BlockLength == 0 ) { return false; }\r
 628 \r
 629     // read in core alignment data, make sure the right size of data was read\r
 630     char x[BAM_CORE_SIZE];\r
 631     if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }\r
 632 \r
 633     if ( IsBigEndian ) {\r
 634         for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) ) { \r
 635           SwapEndian_32p(&x[i]); \r
 636         }\r
 637     }\r
 638     \r
 639     // set BamAlignment 'core' and 'support' data\r
 640     bAlignment.RefID    = BgzfData::UnpackSignedInt(&x[0]);  \r
 641     bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);\r
 642     \r
 643     unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);\r
 644     bAlignment.Bin        = tempValue >> 16;\r
 645     bAlignment.MapQuality = tempValue >> 8 & 0xff;\r
 646     bAlignment.SupportData.QueryNameLength = tempValue & 0xff;\r
 647 \r
 648     tempValue = BgzfData::UnpackUnsignedInt(&x[12]);\r
 649     bAlignment.AlignmentFlag = tempValue >> 16;\r
 650     bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;\r
 651 \r
 652     bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);\r
 653     bAlignment.MateRefID    = BgzfData::UnpackSignedInt(&x[20]);\r
 654     bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);\r
 655     bAlignment.InsertSize   = BgzfData::UnpackSignedInt(&x[28]);\r
 656     \r
 657     // set BamAlignment length\r
 658     bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;\r
 659     \r
 660     // read in character data - make sure proper data size was read\r
 661     bool readCharDataOK = false;\r
 662     const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;\r
 663     char* allCharData = (char*)calloc(sizeof(char), dataLength);\r
 664     \r
 665     if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) { \r
 666       \r
 667         // store 'allCharData' in supportData structure\r
 668         bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);\r
 669         \r
 670         // set success flag\r
 671         readCharDataOK = true;\r
 672         \r
 673         // save CIGAR ops \r
 674         // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly, \r
 675         // even when BamReader::GetNextAlignmentCore() is called \r
 676         const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;\r
 677         uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);\r
 678         CigarOp op;\r
 679         bAlignment.CigarData.clear();\r
 680         bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);\r
 681         for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {\r
 682 \r
 683             // swap if necessary\r
 684             if ( IsBigEndian ) SwapEndian_32(cigarData[i]);\r
 685           \r
 686             // build CigarOp structure\r
 687             op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);\r
 688             op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];\r
 689 \r
 690             // save CigarOp\r
 691             bAlignment.CigarData.push_back(op);\r
 692         }\r
 693     }\r
 694 \r
 695     free(allCharData);\r
 696     return readCharDataOK;\r
 697 }\r
 698 \r
 699 // loads reference data from BAM file\r
 700 void BamReader::BamReaderPrivate::LoadReferenceData(void) {\r
 701 \r
 702     // get number of reference sequences\r
 703     char buffer[4];\r
 704     mBGZF.Read(buffer, 4);\r
 705     unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);\r
 706     if ( IsBigEndian ) { SwapEndian_32(numberRefSeqs); }\r
 707     if (numberRefSeqs == 0) { return; }\r
 708     References.reserve((int)numberRefSeqs);\r
 709 \r
 710     // iterate over all references in header\r
 711     for (unsigned int i = 0; i != numberRefSeqs; ++i) {\r
 712 \r
 713         // get length of reference name\r
 714         mBGZF.Read(buffer, 4);\r
 715         unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);\r
 716         if ( IsBigEndian ) { SwapEndian_32(refNameLength); }\r
 717         char* refName = (char*)calloc(refNameLength, 1);\r
 718 \r
 719         // get reference name and reference sequence length\r
 720         mBGZF.Read(refName, refNameLength);\r
 721         mBGZF.Read(buffer, 4);\r
 722         int refLength = BgzfData::UnpackSignedInt(buffer);\r
 723         if ( IsBigEndian ) { SwapEndian_32(refLength); }\r
 724 \r
 725         // store data for reference\r
 726         RefData aReference;\r
 727         aReference.RefName   = (string)((const char*)refName);\r
 728         aReference.RefLength = refLength;\r
 729         References.push_back(aReference);\r
 730 \r
 731         // clean up calloc-ed temp variable\r
 732         free(refName);\r
 733     }\r
 734 }\r
 735 \r
 736 // opens BAM file (and index)\r
 737 bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename) {\r
 738 \r
 739     Filename = filename;\r
 740     IndexFilename = indexFilename;\r
 741 \r
 742     // open the BGZF file for reading, return false on failure\r
 743     if ( !mBGZF.Open(filename, "rb") ) \r
 744         return false;\r
 745     \r
 746     // retrieve header text & reference data\r
 747     LoadHeaderData();\r
 748     LoadReferenceData();\r
 749 \r
 750     // store file offset of first alignment\r
 751     AlignmentsBeginOffset = mBGZF.Tell();\r
 752 \r
 753     // open index file & load index data (if exists)\r
 754     if ( !IndexFilename.empty() )\r
 755         LoadIndex();\r
 756     \r
 757     // return success\r
 758     return true;\r
 759 }\r
 760 \r
 761 // returns BAM file pointer to beginning of alignment data\r
 762 bool BamReader::BamReaderPrivate::Rewind(void) {\r
 763    \r
 764     // rewind to first alignment\r
 765     if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;\r
 766   \r
 767     // retrieve first alignment data\r
 768     BamAlignment al;\r
 769     if ( !LoadNextAlignment(al) ) return false;\r
 770       \r
 771     // reset default region info using first alignment in file\r
 772     Region.LeftRefID      = al.RefID;\r
 773     Region.LeftPosition   = al.Position;\r
 774     Region.RightRefID     = -1;\r
 775     Region.RightPosition  = -1;\r
 776     IsLeftBoundSpecified  = false;\r
 777     IsRightBoundSpecified = false; \r
 778 \r
 779     // rewind back to before first alignment\r
 780     // return success/fail of seek\r
 781     return mBGZF.Seek(AlignmentsBeginOffset);\r
 782 }\r
 783 \r
 784 // sets a region of interest (with left & right bound reference/position)\r
 785 // attempts a Jump() to left bound as well\r
 786 // returns success/failure of Jump()\r
 787 bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {\r
 788     \r
 789     // save region of interest\r
 790     Region = region;\r
 791     \r
 792     // set flags\r
 793     if ( region.LeftRefID >= 0 && region.LeftPosition >= 0 ) \r
 794         IsLeftBoundSpecified = true;\r
 795     if ( region.RightRefID >= 0 && region.RightPosition >= 0 ) \r
 796         IsRightBoundSpecified = true;\r
 797     \r
 798     // attempt jump to beginning of region, return success/fail of Jump()\r
 799     return Jump( Region.LeftRefID, Region.LeftPosition );\r
 800 }\r