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