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