src/api/internal/BamReader_p.cpp

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