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