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: 10 May 2011 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides the basic functionality for reading BAM files
   9 // ***************************************************************************
  10
  11 #include <api/BamConstants.h>
  12 #include <api/BamReader.h>
  13 #include <api/IBamIODevice.h>
  14 #include <api/internal/BamDeviceFactory_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 <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 void BamReaderPrivate::Close(void) {
  45
  46     // clear header & reference data
  47     m_references.clear();
  48     m_header.Clear();
  49
  50     // close internal
  51     m_randomAccessController.Close();
  52     m_stream.Close();
  53
  54     // clear filename
  55     m_filename.clear();
  56 }
  57
  58 // creates an index file of requested type on current BAM file
  59 bool BamReaderPrivate::CreateIndex(const BamIndex::IndexType& type) {
  60     if ( !IsOpen() ) return false;
  61     return m_randomAccessController.CreateIndex(this, type);
  62 }
  63
  64 // return path & filename of current BAM file
  65 const string BamReaderPrivate::Filename(void) const {
  66     return m_filename;
  67 }
  68
  69 // return header data as std::string
  70 string BamReaderPrivate::GetHeaderText(void) const {
  71     return m_header.ToString();
  72 }
  73
  74 // return header data as SamHeader object
  75 SamHeader BamReaderPrivate::GetSamHeader(void) const {
  76     return m_header.ToSamHeader();
  77 }
  78
  79 // get next alignment (with character data fully parsed)
  80 bool BamReaderPrivate::GetNextAlignment(BamAlignment& alignment) {
  81
  82     // if valid alignment found
  83     if ( GetNextAlignmentCore(alignment) ) {
  84
  85         // store alignment's "source" filename
  86         alignment.Filename = m_filename;
  87
  88         // return success/failure of parsing char data
  89         return alignment.BuildCharData();
  90     }
  91
  92     // no valid alignment found
  93     return false;
  94 }
  95
  96 // retrieves next available alignment core data (returns success/fail)
  97 // ** DOES NOT populate any character data fields (read name, bases, qualities, tag data, filename)
  98 //    these can be accessed, if necessary, from the supportData
  99 // useful for operations requiring ONLY positional or other alignment-related information
 100 bool BamReaderPrivate::GetNextAlignmentCore(BamAlignment& alignment) {
 101
 102     if ( !m_stream.IsOpen() )
 103         return false;
 104
 105     // skip if region is set but has no alignments
 106     if ( m_randomAccessController.HasRegion() &&
 107          !m_randomAccessController.RegionHasAlignments() )
 108     {
 109         return false;
 110     }
 111
 112     // if can't read next alignment
 113     if ( !LoadNextAlignment(alignment) )
 114         return false;
 115
 116     // check alignment's region-overlap state
 117     BamRandomAccessController::RegionState state = m_randomAccessController.AlignmentState(alignment);
 118
 119     // if alignment starts after region, no need to keep reading
 120     if ( state == BamRandomAccessController::AfterRegion )
 121         return false;
 122
 123     // read until overlap is found
 124     while ( state != BamRandomAccessController::OverlapsRegion ) {
 125
 126         // if can't read next alignment
 127         if ( !LoadNextAlignment(alignment) )
 128             return false;
 129
 130         // check alignment's region-overlap state
 131         state = m_randomAccessController.AlignmentState(alignment);
 132
 133         // if alignment starts after region, no need to keep reading
 134         if ( state == BamRandomAccessController::AfterRegion )
 135             return false;
 136     }
 137
 138     // if we get here, we found the next 'valid' alignment
 139     // (e.g. overlaps current region if one was set, simply the next alignment if not)
 140     alignment.SupportData.HasCoreOnly = true;
 141     return true;
 142 }
 143
 144 int BamReaderPrivate::GetReferenceCount(void) const {
 145     return m_references.size();
 146 }
 147
 148 const RefVector& BamReaderPrivate::GetReferenceData(void) const {
 149     return m_references;
 150 }
 151
 152 // returns RefID for given RefName (returns References.size() if not found)
 153 int BamReaderPrivate::GetReferenceID(const string& refName) const {
 154
 155     // retrieve names from reference data
 156     vector<string> refNames;
 157     RefVector::const_iterator refIter = m_references.begin();
 158     RefVector::const_iterator refEnd  = m_references.end();
 159     for ( ; refIter != refEnd; ++refIter)
 160         refNames.push_back( (*refIter).RefName );
 161
 162     // return 'index-of' refName (or -1 if not found)
 163     int index = distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
 164     if ( index == (int)m_references.size() ) return -1;
 165     else return index;
 166 }
 167
 168 bool BamReaderPrivate::HasIndex(void) const {
 169     return m_randomAccessController.HasIndex();
 170 }
 171
 172 bool BamReaderPrivate::IsOpen(void) const {
 173     return m_stream.IsOpen();
 174 }
 175
 176 // load BAM header data
 177 bool BamReaderPrivate::LoadHeaderData(void) {
 178      return m_header.Load(&m_stream);
 179 }
 180
 181 // populates BamAlignment with alignment data under file pointer, returns success/fail
 182 bool BamReaderPrivate::LoadNextAlignment(BamAlignment& alignment) {
 183
 184     // read in the 'block length' value, make sure it's not zero
 185     char buffer[sizeof(uint32_t)];
 186     m_stream.Read(buffer, sizeof(uint32_t));
 187     alignment.SupportData.BlockLength = BamTools::UnpackUnsignedInt(buffer);
 188     if ( m_isBigEndian ) BamTools::SwapEndian_32(alignment.SupportData.BlockLength);
 189     if ( alignment.SupportData.BlockLength == 0 ) return false;
 190
 191     // read in core alignment data, make sure the right size of data was read
 192     char x[Constants::BAM_CORE_SIZE];
 193     if ( m_stream.Read(x, Constants::BAM_CORE_SIZE) != Constants::BAM_CORE_SIZE )
 194         return false;
 195
 196     // swap core endian-ness if necessary
 197     if ( m_isBigEndian ) {
 198         for ( unsigned int i = 0; i < Constants::BAM_CORE_SIZE; i+=sizeof(uint32_t) )
 199             BamTools::SwapEndian_32p(&x[i]);
 200     }
 201
 202     // set BamAlignment 'core' and 'support' data
 203     alignment.RefID    = BamTools::UnpackSignedInt(&x[0]);
 204     alignment.Position = BamTools::UnpackSignedInt(&x[4]);
 205
 206     unsigned int tempValue = BamTools::UnpackUnsignedInt(&x[8]);
 207     alignment.Bin        = tempValue >> 16;
 208     alignment.MapQuality = tempValue >> 8 & 0xff;
 209     alignment.SupportData.QueryNameLength = tempValue & 0xff;
 210
 211     tempValue = BamTools::UnpackUnsignedInt(&x[12]);
 212     alignment.AlignmentFlag = tempValue >> 16;
 213     alignment.SupportData.NumCigarOperations = tempValue & 0xffff;
 214
 215     alignment.SupportData.QuerySequenceLength = BamTools::UnpackUnsignedInt(&x[16]);
 216     alignment.MateRefID    = BamTools::UnpackSignedInt(&x[20]);
 217     alignment.MatePosition = BamTools::UnpackSignedInt(&x[24]);
 218     alignment.InsertSize   = BamTools::UnpackSignedInt(&x[28]);
 219
 220     // set BamAlignment length
 221     alignment.Length = alignment.SupportData.QuerySequenceLength;
 222
 223     // read in character data - make sure proper data size was read
 224     bool readCharDataOK = false;
 225     const unsigned int dataLength = alignment.SupportData.BlockLength - Constants::BAM_CORE_SIZE;
 226     char* allCharData = (char*)calloc(sizeof(char), dataLength);
 227
 228     if ( m_stream.Read(allCharData, dataLength) == dataLength ) {
 229
 230         // store 'allCharData' in supportData structure
 231         alignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
 232
 233         // set success flag
 234         readCharDataOK = true;
 235
 236         // save CIGAR ops
 237         // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly,
 238         // even when GetNextAlignmentCore() is called
 239         const unsigned int cigarDataOffset = alignment.SupportData.QueryNameLength;
 240         uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
 241         CigarOp op;
 242         alignment.CigarData.clear();
 243         alignment.CigarData.reserve(alignment.SupportData.NumCigarOperations);
 244         for ( unsigned int i = 0; i < alignment.SupportData.NumCigarOperations; ++i ) {
 245
 246             // swap endian-ness if necessary
 247             if ( m_isBigEndian ) BamTools::SwapEndian_32(cigarData[i]);
 248
 249             // build CigarOp structure
 250             op.Length = (cigarData[i] >> Constants::BAM_CIGAR_SHIFT);
 251             op.Type   = Constants::BAM_CIGAR_LOOKUP[ (cigarData[i] & Constants::BAM_CIGAR_MASK) ];
 252
 253             // save CigarOp
 254             alignment.CigarData.push_back(op);
 255         }
 256     }
 257
 258     // clean up & return parsing success/failure
 259     free(allCharData);
 260     return readCharDataOK;
 261 }
 262
 263 // loads reference data from BAM file
 264 bool BamReaderPrivate::LoadReferenceData(void) {
 265
 266     // get number of reference sequences
 267     char buffer[sizeof(uint32_t)];
 268     m_stream.Read(buffer, sizeof(uint32_t));
 269     uint32_t numberRefSeqs = BamTools::UnpackUnsignedInt(buffer);
 270     if ( m_isBigEndian ) BamTools::SwapEndian_32(numberRefSeqs);
 271     m_references.reserve((int)numberRefSeqs);
 272
 273     // iterate over all references in header
 274     for ( unsigned int i = 0; i != numberRefSeqs; ++i ) {
 275
 276         // get length of reference name
 277         m_stream.Read(buffer, sizeof(uint32_t));
 278         uint32_t refNameLength = BamTools::UnpackUnsignedInt(buffer);
 279         if ( m_isBigEndian ) BamTools::SwapEndian_32(refNameLength);
 280         char* refName = (char*)calloc(refNameLength, 1);
 281
 282         // get reference name and reference sequence length
 283         m_stream.Read(refName, refNameLength);
 284         m_stream.Read(buffer, sizeof(int32_t));
 285         int32_t refLength = BamTools::UnpackSignedInt(buffer);
 286         if ( m_isBigEndian ) BamTools::SwapEndian_32(refLength);
 287
 288         // store data for reference
 289         RefData aReference;
 290         aReference.RefName   = (string)((const char*)refName);
 291         aReference.RefLength = refLength;
 292         m_references.push_back(aReference);
 293
 294         // clean up calloc-ed temp variable
 295         free(refName);
 296     }
 297
 298     // return success
 299     return true;
 300 }
 301
 302 bool BamReaderPrivate::LocateIndex(const BamIndex::IndexType& preferredType) {
 303     return m_randomAccessController.LocateIndex(this, preferredType);
 304 }
 305
 306 // opens BAM file (and index)
 307 bool BamReaderPrivate::Open(const string& filename) {
 308
 309     // make sure we're starting with a fresh slate
 310     Close();
 311
 312     // attempt to open BgzfStream for reading
 313     if ( !m_stream.Open(filename, IBamIODevice::ReadOnly) ) {
 314         cerr << "BamReader ERROR: Could not open BGZF stream for " << filename << endl;
 315         return false;
 316     }
 317
 318     // attempt to load header data
 319     if ( !LoadHeaderData() ) {
 320         cerr << "BamReader ERROR: Could not load header data for " << filename << endl;
 321         Close();
 322         return false;
 323     }
 324
 325     // attempt to load reference data
 326     if ( !LoadReferenceData() ) {
 327         cerr << "BamReader ERROR: Could not load reference data for " << filename << endl;
 328         Close();
 329         return false;
 330     }
 331
 332     // if all OK, store filename & offset of first alignment
 333     m_filename = filename;
 334     m_alignmentsBeginOffset = m_stream.Tell();
 335
 336     // return success
 337     return true;
 338 }
 339
 340 bool BamReaderPrivate::OpenIndex(const std::string& indexFilename) {
 341     return m_randomAccessController.OpenIndex(indexFilename, this);
 342 }
 343
 344 // returns BAM file pointer to beginning of alignment data
 345 bool BamReaderPrivate::Rewind(void) {
 346
 347     if ( !m_stream.IsOpen() ) {
 348         cerr << "BRP::Rewind() - stream not open" << endl;
 349         return false;
 350     }
 351
 352     // attempt rewind to first alignment
 353     if ( !m_stream.Seek(m_alignmentsBeginOffset) ) {
 354         cerr << "BRP::Rewind() - could not seek to ABO (1st time)" << endl;
 355         return false;
 356     }
 357
 358     // verify that we can read first alignment
 359     BamAlignment al;
 360     if ( !LoadNextAlignment(al) ) {
 361         cerr << "BRP::Rewind() - could not read first alignment" << endl;
 362         return false;
 363     }
 364
 365     // reset region
 366     m_randomAccessController.ClearRegion();
 367
 368     // rewind back to beginning of first alignment
 369     // return success/fail of seek
 370     return m_stream.Seek(m_alignmentsBeginOffset);
 371 }
 372
 373 bool BamReaderPrivate::Seek(const int64_t& position) {
 374     return m_stream.Seek(position);
 375 }
 376
 377 void BamReaderPrivate::SetIndex(BamIndex* index) {
 378     m_randomAccessController.SetIndex(index);
 379 }
 380
 381 // change the index caching behavior
 382 void BamReaderPrivate::SetIndexCacheMode(const BamIndex::IndexCacheMode& mode) {
 383     m_randomAccessController.SetIndexCacheMode(mode);
 384 }
 385
 386 // sets current region & attempts to jump to it
 387 // returns success/failure
 388 bool BamReaderPrivate::SetRegion(const BamRegion& region) {
 389     return m_randomAccessController.SetRegion(this, region, m_references.size());
 390 }
 391
 392 int64_t BamReaderPrivate::Tell(void) const {
 393     return m_stream.Tell();
 394 }