src/api/BamReader.cpp

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