src/api/BamReader.cpp

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