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