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