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