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