src/api/BamReader.cpp

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