src/api/internal/BamReader_p.cpp

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