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 <<<<<<< HEAD:src/api/internal/BamReader_p.cpp
   6 // Last modified: 14 November 2011 (DB)
   7 =======
   8 // Last modified: 25 October 2011 (DB)
   9 >>>>>>> remoteio:src/api/internal/bam/BamReader_p.cpp
  10 // ---------------------------------------------------------------------------
  11 // Provides the basic functionality for reading BAM files
  12 // ***************************************************************************
  13
  14 #include "api/BamConstants.h"
  15 #include "api/BamReader.h"
  16 #include "api/IBamIODevice.h"
  17 #include "api/internal/bam/BamHeader_p.h"
  18 #include "api/internal/bam/BamRandomAccessController_p.h"
  19 #include "api/internal/bam/BamReader_p.h"
  20 #include "api/internal/index/BamStandardIndex_p.h"
  21 #include "api/internal/index/BamToolsIndex_p.h"
  22 #include "api/internal/io/BamDeviceFactory_p.h"
  23 #include "api/internal/utils/BamException_p.h"
  24 using namespace BamTools;
  25 using namespace BamTools::Internal;
  26
  27 #include <algorithm>
  28 #include <cassert>
  29 #include <iostream>
  30 #include <iterator>
  31 #include <vector>
  32 using namespace std;
  33
  34 // constructor
  35 BamReaderPrivate::BamReaderPrivate(BamReader* parent)
  36     : m_alignmentsBeginOffset(0)
  37     , m_parent(parent)
  38 {
  39     m_isBigEndian = BamTools::SystemIsBigEndian();
  40 }
  41
  42 // destructor
  43 BamReaderPrivate::~BamReaderPrivate(void) {
  44     Close();
  45 }
  46
  47 // closes the BAM file
  48 bool BamReaderPrivate::Close(void) {
  49
  50     // clear BAM metadata
  51     m_references.clear();
  52     m_header.Clear();
  53
  54     // clear filename
  55     m_filename.clear();
  56
  57     // close random access controller
  58     m_randomAccessController.Close();
  59
  60     // if stream is open, attempt close
  61     if ( IsOpen() ) {
  62         try {
  63             m_stream.Close();
  64         } catch ( BamException& e ) {
  65             const string streamError = e.what();
  66             const string message = string("encountered error closing BAM file: \n\t") + streamError;
  67             SetErrorString("BamReader::Close", message);
  68             return false;
  69         }
  70     }
  71
  72     // return success
  73     return true;
  74 }
  75
  76 // creates an index file of requested type on current BAM file
  77 bool BamReaderPrivate::CreateIndex(const BamIndex::IndexType& type) {
  78
  79     // skip if BAM file not open
  80     if ( !IsOpen() ) {
  81         SetErrorString("BamReader::CreateIndex", "cannot create index on unopened BAM file");
  82         return false;
  83     }
  84
  85     // attempt to create index
  86     if ( m_randomAccessController.CreateIndex(this, type) )
  87         return true;
  88     else {
  89         const string bracError = m_randomAccessController.GetErrorString();
  90         const string message = string("could not create index: \n\t") + bracError;
  91         SetErrorString("BamReader::CreateIndex", message);
  92         return false;
  93     }
  94 }
  95
  96 // return path & filename of current BAM file
  97 const string BamReaderPrivate::Filename(void) const {
  98     return m_filename;
  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     m_stream.Read(buffer, sizeof(uint32_t));
 240     alignment.SupportData.BlockLength = BamTools::UnpackUnsignedInt(buffer);
 241     if ( m_isBigEndian ) BamTools::SwapEndian_32(alignment.SupportData.BlockLength);
 242     if ( alignment.SupportData.BlockLength == 0 )
 243         return false;
 244
 245     // read in core alignment data, make sure the right size of data was read
 246     char x[Constants::BAM_CORE_SIZE];
 247     if ( m_stream.Read(x, Constants::BAM_CORE_SIZE) != Constants::BAM_CORE_SIZE )
 248         return false;
 249
 250     // swap core endian-ness if necessary
 251     if ( m_isBigEndian ) {
 252         for ( unsigned int i = 0; i < Constants::BAM_CORE_SIZE; i+=sizeof(uint32_t) )
 253             BamTools::SwapEndian_32p(&x[i]);
 254     }
 255
 256     // set BamAlignment 'core' and 'support' data
 257     alignment.RefID    = BamTools::UnpackSignedInt(&x[0]);
 258     alignment.Position = BamTools::UnpackSignedInt(&x[4]);
 259
 260     unsigned int tempValue = BamTools::UnpackUnsignedInt(&x[8]);
 261     alignment.Bin        = tempValue >> 16;
 262     alignment.MapQuality = tempValue >> 8 & 0xff;
 263     alignment.SupportData.QueryNameLength = tempValue & 0xff;
 264
 265     tempValue = BamTools::UnpackUnsignedInt(&x[12]);
 266     alignment.AlignmentFlag = tempValue >> 16;
 267     alignment.SupportData.NumCigarOperations = tempValue & 0xffff;
 268
 269     alignment.SupportData.QuerySequenceLength = BamTools::UnpackUnsignedInt(&x[16]);
 270     alignment.MateRefID    = BamTools::UnpackSignedInt(&x[20]);
 271     alignment.MatePosition = BamTools::UnpackSignedInt(&x[24]);
 272     alignment.InsertSize   = BamTools::UnpackSignedInt(&x[28]);
 273
 274     // set BamAlignment length
 275     alignment.Length = alignment.SupportData.QuerySequenceLength;
 276
 277     // read in character data - make sure proper data size was read
 278     bool readCharDataOK = false;
 279     const unsigned int dataLength = alignment.SupportData.BlockLength - Constants::BAM_CORE_SIZE;
 280     RaiiBuffer allCharData(dataLength);
 281
 282     if ( m_stream.Read(allCharData.Buffer, dataLength) == dataLength ) {
 283
 284         // store 'allCharData' in supportData structure
 285         alignment.SupportData.AllCharData.assign((const char*)allCharData.Buffer, dataLength);
 286
 287         // set success flag
 288         readCharDataOK = true;
 289
 290         // save CIGAR ops
 291         // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly,
 292         // even when GetNextAlignmentCore() is called
 293         const unsigned int cigarDataOffset = alignment.SupportData.QueryNameLength;
 294         uint32_t* cigarData = (uint32_t*)(allCharData.Buffer + cigarDataOffset);
 295         CigarOp op;
 296         alignment.CigarData.clear();
 297         alignment.CigarData.reserve(alignment.SupportData.NumCigarOperations);
 298         for ( unsigned int i = 0; i < alignment.SupportData.NumCigarOperations; ++i ) {
 299
 300             // swap endian-ness if necessary
 301             if ( m_isBigEndian ) BamTools::SwapEndian_32(cigarData[i]);
 302
 303             // build CigarOp structure
 304             op.Length = (cigarData[i] >> Constants::BAM_CIGAR_SHIFT);
 305             op.Type   = Constants::BAM_CIGAR_LOOKUP[ (cigarData[i] & Constants::BAM_CIGAR_MASK) ];
 306
 307             // save CigarOp
 308             alignment.CigarData.push_back(op);
 309         }
 310     }
 311
 312     // return success/failure
 313     return readCharDataOK;
 314 }
 315
 316 // loads reference data from BAM file
 317 bool BamReaderPrivate::LoadReferenceData(void) {
 318
 319     // get number of reference sequences
 320     char buffer[sizeof(uint32_t)];
 321     m_stream.Read(buffer, sizeof(uint32_t));
 322     uint32_t numberRefSeqs = BamTools::UnpackUnsignedInt(buffer);
 323     if ( m_isBigEndian ) BamTools::SwapEndian_32(numberRefSeqs);
 324     m_references.reserve((int)numberRefSeqs);
 325
 326     // iterate over all references in header
 327     for ( unsigned int i = 0; i != numberRefSeqs; ++i ) {
 328
 329         // get length of reference name
 330         m_stream.Read(buffer, sizeof(uint32_t));
 331         uint32_t refNameLength = BamTools::UnpackUnsignedInt(buffer);
 332         if ( m_isBigEndian ) BamTools::SwapEndian_32(refNameLength);
 333         RaiiBuffer refName(refNameLength);
 334
 335         // get reference name and reference sequence length
 336         m_stream.Read(refName.Buffer, refNameLength);
 337         m_stream.Read(buffer, sizeof(int32_t));
 338         int32_t refLength = BamTools::UnpackSignedInt(buffer);
 339         if ( m_isBigEndian ) BamTools::SwapEndian_32(refLength);
 340
 341         // store data for reference
 342         RefData aReference;
 343         aReference.RefName   = (string)((const char*)refName.Buffer);
 344         aReference.RefLength = refLength;
 345         m_references.push_back(aReference);
 346     }
 347
 348     // return success
 349     return true;
 350 }
 351
 352 bool BamReaderPrivate::LocateIndex(const BamIndex::IndexType& preferredType) {
 353
 354     if ( m_randomAccessController.LocateIndex(this, preferredType) )
 355         return true;
 356     else {
 357         const string bracError = m_randomAccessController.GetErrorString();
 358         const string message = string("could not locate index: \n\t") + bracError;
 359         SetErrorString("BamReader::LocateIndex", message);
 360         return false;
 361     }
 362 }
 363
 364 // opens BAM file (and index)
 365 bool BamReaderPrivate::Open(const string& filename) {
 366
 367     try {
 368
 369         // make sure we're starting with fresh state
 370         Close();
 371
 372         // open BgzfStream
 373         m_stream.Open(filename, IBamIODevice::ReadOnly);
 374
 375         // load BAM metadata
 376         LoadHeaderData();
 377         LoadReferenceData();
 378
 379         // store filename & offset of first alignment
 380         m_filename = filename;
 381         m_alignmentsBeginOffset = m_stream.Tell();
 382
 383         // return success
 384         return true;
 385
 386     } catch ( BamException& e ) {
 387         const string error = e.what();
 388         const string message = string("could not open file: ") + filename +
 389                                "\n\t" + error;
 390         SetErrorString("BamReader::Open", message);
 391         return false;
 392     }
 393 }
 394
 395 bool BamReaderPrivate::OpenIndex(const std::string& indexFilename) {
 396
 397     if ( m_randomAccessController.OpenIndex(indexFilename, this) )
 398         return true;
 399     else {
 400         const string bracError = m_randomAccessController.GetErrorString();
 401         const string message = string("could not open index: \n\t") + bracError;
 402         SetErrorString("BamReader::OpenIndex", message);
 403         return false;
 404     }
 405 }
 406
 407 // returns BAM file pointer to beginning of alignment data
 408 bool BamReaderPrivate::Rewind(void) {
 409
 410     // reset region
 411     m_randomAccessController.ClearRegion();
 412
 413     // return status of seeking back to first alignment
 414     if ( Seek(m_alignmentsBeginOffset) )
 415         return true;
 416     else {
 417         const string currentError = m_errorString;
 418         const string message = string("could not rewind: \n\t") + currentError;
 419         SetErrorString("BamReader::Rewind", message);
 420         return false;
 421     }
 422 }
 423
 424 bool BamReaderPrivate::Seek(const int64_t& position) {
 425
 426     // skip if BAM file not open
 427     if ( !IsOpen() ) {
 428         SetErrorString("BamReader::Seek", "cannot seek on unopened BAM file");
 429         return false;
 430     }
 431
 432     try {
 433         m_stream.Seek(position);
 434         return true;
 435     }
 436     catch ( BamException& e ) {
 437         const string streamError = e.what();
 438         const string message = string("could not seek in BAM file: \n\t") + streamError;
 439         SetErrorString("BamReader::Seek", message);
 440         return false;
 441     }
 442 }
 443
 444 void BamReaderPrivate::SetErrorString(const string& where, const string& what) {
 445     static const string SEPARATOR = ": ";
 446     m_errorString = where + SEPARATOR + what;
 447 }
 448
 449 void BamReaderPrivate::SetIndex(BamIndex* index) {
 450     m_randomAccessController.SetIndex(index);
 451 }
 452
 453 // sets current region & attempts to jump to it
 454 // returns success/failure
 455 bool BamReaderPrivate::SetRegion(const BamRegion& region) {
 456
 457     if ( m_randomAccessController.SetRegion(region, m_references.size()) )
 458         return true;
 459     else {
 460         const string bracError = m_randomAccessController.GetErrorString();
 461         const string message = string("could not set region: \n\t") + bracError;
 462         SetErrorString("BamReader::SetRegion", message);
 463         return false;
 464     }
 465 }
 466
 467 int64_t BamReaderPrivate::Tell(void) const {
 468     return m_stream.Tell();
 469 }