src/api/internal/bam/BamReader_p.cpp

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