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