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