src/api/internal/BamReader_p.cpp

   1 // ***************************************************************************
   2 // BamReader_p.cpp (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 7 October 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides the basic functionality for reading BAM files
   8 // ***************************************************************************
   9
  10 #include <api/BamConstants.h>
  11 #include <api/BamReader.h>
  12 #include <api/internal/BamException_p.h>
  13 #include <api/internal/BamHeader_p.h>
  14 #include <api/internal/BamRandomAccessController_p.h>
  15 #include <api/internal/BamReader_p.h>
  16 #include <api/internal/BamStandardIndex_p.h>
  17 #include <api/internal/BamToolsIndex_p.h>
  18 #include <api/internal/BgzfStream_p.h>
  19 using namespace BamTools;
  20 using namespace BamTools::Internal;
  21
  22 #include <algorithm>
  23 #include <cassert>
  24 #include <iostream>
  25 #include <iterator>
  26 #include <vector>
  27 using namespace std;
  28
  29 // constructor
  30 BamReaderPrivate::BamReaderPrivate(BamReader* parent)
  31     : m_alignmentsBeginOffset(0)
  32     , m_parent(parent)
  33 {
  34     m_isBigEndian = BamTools::SystemIsBigEndian();
  35 }
  36
  37 // destructor
  38 BamReaderPrivate::~BamReaderPrivate(void) {
  39     Close();
  40 }
  41
  42 // closes the BAM file
  43 bool BamReaderPrivate::Close(void) {
  44
  45     // clear BAM metadata
  46     m_references.clear();
  47     m_header.Clear();
  48
  49     // clear filename
  50     m_filename.clear();
  51
  52     // close random access controller
  53     m_randomAccessController.Close();
  54
  55     // if stream is open, attempt close
  56     if ( IsOpen() ) {
  57         try {
  58             m_stream.Close();
  59         } catch ( BamException& e ) {
  60             const string streamError = e.what();
  61             const string message = string("encountered error closing BAM file: \n\t") + streamError;
  62             SetErrorString("BamReader::Close", message);
  63             return false;
  64         }
  65     }
  66
  67     // return success
  68     return true;
  69 }
  70
  71 // creates an index file of requested type on current BAM file
  72 bool BamReaderPrivate::CreateIndex(const BamIndex::IndexType& type) {
  73
  74     // skip if BAM file not open
  75     if ( !IsOpen() ) {
  76         SetErrorString("BamReader::CreateIndex", "cannot create index on unopened BAM file");
  77         return false;
  78     }
  79
  80     // attempt to create index
  81     if ( m_randomAccessController.CreateIndex(this, type) )
  82         return true;
  83     else {
  84         const string bracError = m_randomAccessController.GetErrorString();
  85         const string message = string("could not create index: \n\t") + bracError;
  86         SetErrorString("BamReader::CreateIndex", message);
  87         return false;
  88     }
  89 }
  90
  91 // return path & filename of current BAM file
  92 const string BamReaderPrivate::Filename(void) const {
  93     return m_filename;
  94 }
  95
  96 string BamReaderPrivate::GetErrorString(void) const {
  97     return m_errorString;
  98 }
  99
 100 // return header data as std::string
 101 string BamReaderPrivate::GetHeaderText(void) const {
 102     return m_header.ToString();
 103 }
 104
 105 // return header data as SamHeader object
 106 SamHeader BamReaderPrivate::GetSamHeader(void) const {
 107     return m_header.ToSamHeader();
 108 }
 109
 110 // get next alignment (with character data fully parsed)
 111 bool BamReaderPrivate::GetNextAlignment(BamAlignment& alignment) {
 112
 113     // if valid alignment found
 114     if ( GetNextAlignmentCore(alignment) ) {
 115
 116         // store alignment's "source" filename
 117         alignment.Filename = m_filename;
 118
 119         // return success/failure of parsing char data
 120         if ( alignment.BuildCharData() )
 121             return true;
 122         else {
 123             const string alError = alignment.GetErrorString();
 124             const string message = string("could not populate alignment data: \n\t") + alError;
 125             SetErrorString("BamReader::GetNextAlignment", message);
 126             return false;
 127         }
 128     }
 129
 130     // no valid alignment found
 131     return false;
 132 }
 133
 134 // retrieves next available alignment core data (returns success/fail)
 135 // ** DOES NOT populate any character data fields (read name, bases, qualities, tag data, filename)
 136 //    these can be accessed, if necessary, from the supportData
 137 // useful for operations requiring ONLY positional or other alignment-related information
 138 bool BamReaderPrivate::GetNextAlignmentCore(BamAlignment& alignment) {
 139
 140     try {
 141
 142         // skip if region is set but has no alignments
 143         if ( m_randomAccessController.HasRegion() &&
 144              !m_randomAccessController.RegionHasAlignments() )
 145         {
 146             return false;
 147         }
 148
 149         // if can't read next alignment
 150         if ( !LoadNextAlignment(alignment) )
 151             return false;
 152
 153         // check alignment's region-overlap state
 154         BamRandomAccessController::RegionState state = m_randomAccessController.AlignmentState(alignment);
 155
 156         // if alignment starts after region, no need to keep reading
 157         if ( state == BamRandomAccessController::AfterRegion )
 158             return false;
 159
 160         // read until overlap is found
 161         while ( state != BamRandomAccessController::OverlapsRegion ) {
 162
 163             // if can't read next alignment
 164             if ( !LoadNextAlignment(alignment) )
 165                 return false;
 166
 167             // check alignment's region-overlap state
 168             state = m_randomAccessController.AlignmentState(alignment);
 169
 170             // if alignment starts after region, no need to keep reading
 171             if ( state == BamRandomAccessController::AfterRegion )
 172                 return false;
 173         }
 174
 175         // if we get here, we found the next 'valid' alignment
 176         // (e.g. overlaps current region if one was set, simply the next alignment if not)
 177         alignment.SupportData.HasCoreOnly = true;
 178         return true;
 179
 180     } catch ( BamException& e ) {
 181         const string streamError = e.what();
 182         const string message = string("encountered error reading BAM alignment: \n\t") + streamError;
 183         SetErrorString("BamReader::GetNextAlignmentCore", message);
 184         return false;
 185     }
 186 }
 187
 188 int BamReaderPrivate::GetReferenceCount(void) const {
 189     return m_references.size();
 190 }
 191
 192 const RefVector& BamReaderPrivate::GetReferenceData(void) const {
 193     return m_references;
 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 = m_references.begin();
 202     RefVector::const_iterator refEnd  = m_references.end();
 203     for ( ; refIter != refEnd; ++refIter)
 204         refNames.push_back( (*refIter).RefName );
 205
 206     // return 'index-of' refName (or -1 if not found)
 207     int index = distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
 208     if ( index == (int)m_references.size() ) return -1;
 209     else return index;
 210 }
 211
 212 bool BamReaderPrivate::HasIndex(void) const {
 213     return m_randomAccessController.HasIndex();
 214 }
 215
 216 bool BamReaderPrivate::IsOpen(void) const {
 217     return m_stream.IsOpen;
 218 }
 219
 220 // load BAM header data
 221 void BamReaderPrivate::LoadHeaderData(void) {
 222     m_header.Load(&m_stream);
 223 }
 224
 225 // populates BamAlignment with alignment data under file pointer, returns success/fail
 226 bool BamReaderPrivate::LoadNextAlignment(BamAlignment& alignment) {
 227
 228     // read in the 'block length' value, make sure it's not zero
 229     char buffer[sizeof(uint32_t)];
 230     m_stream.Read(buffer, sizeof(uint32_t));
 231     alignment.SupportData.BlockLength = BamTools::UnpackUnsignedInt(buffer);
 232     if ( m_isBigEndian ) BamTools::SwapEndian_32(alignment.SupportData.BlockLength);
 233     if ( alignment.SupportData.BlockLength == 0 )
 234         return false;
 235
 236     // read in core alignment data, make sure the right size of data was read
 237     char x[Constants::BAM_CORE_SIZE];
 238     if ( m_stream.Read(x, Constants::BAM_CORE_SIZE) != Constants::BAM_CORE_SIZE )
 239         return false;
 240
 241     // swap core endian-ness if necessary
 242     if ( m_isBigEndian ) {
 243         for ( int i = 0; i < Constants::BAM_CORE_SIZE; i+=sizeof(uint32_t) )
 244             BamTools::SwapEndian_32p(&x[i]);
 245     }
 246
 247     // set BamAlignment 'core' and 'support' data
 248     alignment.RefID    = BamTools::UnpackSignedInt(&x[0]);
 249     alignment.Position = BamTools::UnpackSignedInt(&x[4]);
 250
 251     unsigned int tempValue = BamTools::UnpackUnsignedInt(&x[8]);
 252     alignment.Bin        = tempValue >> 16;
 253     alignment.MapQuality = tempValue >> 8 & 0xff;
 254     alignment.SupportData.QueryNameLength = tempValue & 0xff;
 255
 256     tempValue = BamTools::UnpackUnsignedInt(&x[12]);
 257     alignment.AlignmentFlag = tempValue >> 16;
 258     alignment.SupportData.NumCigarOperations = tempValue & 0xffff;
 259
 260     alignment.SupportData.QuerySequenceLength = BamTools::UnpackUnsignedInt(&x[16]);
 261     alignment.MateRefID    = BamTools::UnpackSignedInt(&x[20]);
 262     alignment.MatePosition = BamTools::UnpackSignedInt(&x[24]);
 263     alignment.InsertSize   = BamTools::UnpackSignedInt(&x[28]);
 264
 265     // set BamAlignment length
 266     alignment.Length = alignment.SupportData.QuerySequenceLength;
 267
 268     // read in character data - make sure proper data size was read
 269     bool readCharDataOK = false;
 270     const unsigned int dataLength = alignment.SupportData.BlockLength - Constants::BAM_CORE_SIZE;
 271     RaiiBuffer allCharData(dataLength);
 272
 273     if ( m_stream.Read(allCharData.Buffer, dataLength) == dataLength ) {
 274
 275         // store 'allCharData' in supportData structure
 276         alignment.SupportData.AllCharData.assign((const char*)allCharData.Buffer, dataLength);
 277
 278         // set success flag
 279         readCharDataOK = true;
 280
 281         // save CIGAR ops
 282         // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly,
 283         // even when GetNextAlignmentCore() is called
 284         const unsigned int cigarDataOffset = alignment.SupportData.QueryNameLength;
 285         uint32_t* cigarData = (uint32_t*)(allCharData.Buffer + cigarDataOffset);
 286         CigarOp op;
 287         alignment.CigarData.clear();
 288         alignment.CigarData.reserve(alignment.SupportData.NumCigarOperations);
 289         for ( unsigned int i = 0; i < alignment.SupportData.NumCigarOperations; ++i ) {
 290
 291             // swap endian-ness if necessary
 292             if ( m_isBigEndian ) BamTools::SwapEndian_32(cigarData[i]);
 293
 294             // build CigarOp structure
 295             op.Length = (cigarData[i] >> Constants::BAM_CIGAR_SHIFT);
 296             op.Type   = Constants::BAM_CIGAR_LOOKUP[ (cigarData[i] & Constants::BAM_CIGAR_MASK) ];
 297
 298             // save CigarOp
 299             alignment.CigarData.push_back(op);
 300         }
 301     }
 302
 303     // return success/failure
 304     return readCharDataOK;
 305 }
 306
 307 // loads reference data from BAM file
 308 bool BamReaderPrivate::LoadReferenceData(void) {
 309
 310     // get number of reference sequences
 311     char buffer[sizeof(uint32_t)];
 312     m_stream.Read(buffer, sizeof(uint32_t));
 313     uint32_t numberRefSeqs = BamTools::UnpackUnsignedInt(buffer);
 314     if ( m_isBigEndian ) BamTools::SwapEndian_32(numberRefSeqs);
 315     m_references.reserve((int)numberRefSeqs);
 316
 317     // iterate over all references in header
 318     for ( unsigned int i = 0; i != numberRefSeqs; ++i ) {
 319
 320         // get length of reference name
 321         m_stream.Read(buffer, sizeof(uint32_t));
 322         uint32_t refNameLength = BamTools::UnpackUnsignedInt(buffer);
 323         if ( m_isBigEndian ) BamTools::SwapEndian_32(refNameLength);
 324         RaiiBuffer refName(refNameLength);
 325
 326         // get reference name and reference sequence length
 327         m_stream.Read(refName.Buffer, refNameLength);
 328         m_stream.Read(buffer, sizeof(int32_t));
 329         int32_t refLength = BamTools::UnpackSignedInt(buffer);
 330         if ( m_isBigEndian ) BamTools::SwapEndian_32(refLength);
 331
 332         // store data for reference
 333         RefData aReference;
 334         aReference.RefName   = (string)((const char*)refName.Buffer);
 335         aReference.RefLength = refLength;
 336         m_references.push_back(aReference);
 337     }
 338
 339     // return success
 340     return true;
 341 }
 342
 343 bool BamReaderPrivate::LocateIndex(const BamIndex::IndexType& preferredType) {
 344
 345     if ( m_randomAccessController.LocateIndex(this, preferredType) )
 346         return true;
 347     else {
 348         const string bracError = m_randomAccessController.GetErrorString();
 349         const string message = string("could not locate index: \n\t") + bracError;
 350         SetErrorString("BamReader::LocateIndex", message);
 351         return false;
 352     }
 353 }
 354
 355 // opens BAM file (and index)
 356 bool BamReaderPrivate::Open(const string& filename) {
 357
 358     bool result;
 359
 360     try {
 361
 362         // make sure we're starting with fresh state
 363         Close();
 364
 365         // open BgzfStream
 366         m_stream.Open(filename, "rb");
 367         assert(m_stream);
 368
 369         // load BAM metadata
 370         LoadHeaderData();
 371         LoadReferenceData();
 372
 373         // store filename & offset of first alignment
 374         m_filename = filename;
 375         m_alignmentsBeginOffset = m_stream.Tell();
 376
 377         // set flag
 378         result = true;
 379
 380     } catch ( BamException& e ) {
 381         const string error = e.what();
 382         const string message = string("could not open file: ") + filename +
 383                                "\n\t" + error;
 384         SetErrorString("BamReader::Open", message);
 385         return false;
 386     }
 387
 388     // return success/failure
 389     return result;
 390 }
 391
 392 bool BamReaderPrivate::OpenIndex(const std::string& indexFilename) {
 393
 394     if ( m_randomAccessController.OpenIndex(indexFilename, this) )
 395         return true;
 396     else {
 397         const string bracError = m_randomAccessController.GetErrorString();
 398         const string message = string("could not open index: \n\t") + bracError;
 399         SetErrorString("BamReader::OpenIndex", message);
 400         return false;
 401     }
 402 }
 403
 404 // returns BAM file pointer to beginning of alignment data
 405 bool BamReaderPrivate::Rewind(void) {
 406
 407     // reset region
 408     m_randomAccessController.ClearRegion();
 409
 410     // return status of seeking back to first alignment
 411     if ( Seek(m_alignmentsBeginOffset) )
 412         return true;
 413     else {
 414         const string currentError = m_errorString;
 415         const string message = string("could not rewind: \n\t") + currentError;
 416         SetErrorString("BamReader::Rewind", message);
 417         return false;
 418     }
 419 }
 420
 421 bool BamReaderPrivate::Seek(const int64_t& position) {
 422
 423     // skip if BAM file not open
 424     if ( !IsOpen() ) {
 425         SetErrorString("BamReader::Seek", "cannot seek on unopened BAM file");
 426         return false;
 427     }
 428
 429     try {
 430         m_stream.Seek(position);
 431         return true;
 432     }
 433     catch ( BamException& e ) {
 434         const string streamError = e.what();
 435         const string message = string("could not seek in BAM file: \n\t") + streamError;
 436         SetErrorString("BamReader::Seek", message);
 437         return false;
 438     }
 439 }
 440
 441 void BamReaderPrivate::SetErrorString(const string& where, const string& what) {
 442     static const string SEPARATOR = ": ";
 443     m_errorString = where + SEPARATOR + what;
 444 }
 445
 446 void BamReaderPrivate::SetIndex(BamIndex* index) {
 447     m_randomAccessController.SetIndex(index);
 448 }
 449
 450 // change the index caching behavior
 451 void BamReaderPrivate::SetIndexCacheMode(const BamIndex::IndexCacheMode& mode) {
 452     m_randomAccessController.SetIndexCacheMode(mode);
 453 }
 454
 455 // sets current region & attempts to jump to it
 456 // returns success/failure
 457 bool BamReaderPrivate::SetRegion(const BamRegion& region) {
 458
 459     if ( m_randomAccessController.SetRegion(region, m_references.size()) )
 460         return true;
 461     else {
 462         const string bracError = m_randomAccessController.GetErrorString();
 463         const string message = string("could not set region: \n\t") + bracError;
 464         SetErrorString("BamReader::SetRegion", message);
 465         return false;
 466     }
 467 }
 468
 469 int64_t BamReaderPrivate::Tell(void) const {
 470     return m_stream.Tell();
 471 }