X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fapi%2Finternal%2FBamStandardIndex_p.cpp;h=c492899a20167d4924be71ff25d18c0b55fc0a18;hb=9f1ce8c47aeadb6dc1320b52ee671c3341b97935;hp=cf0e2c1b0a2006d3ab33ebe5414fec45efab0a41;hpb=8c80d760637f8df39262683cd2570f0589423d36;p=bamtools.git diff --git a/src/api/internal/BamStandardIndex_p.cpp b/src/api/internal/BamStandardIndex_p.cpp index cf0e2c1..c492899 100644 --- a/src/api/internal/BamStandardIndex_p.cpp +++ b/src/api/internal/BamStandardIndex_p.cpp @@ -1,18 +1,16 @@ // *************************************************************************** // BamStandardIndex.cpp (c) 2010 Derek Barnett // Marth Lab, Department of Biology, Boston College -// All rights reserved. // --------------------------------------------------------------------------- -// Last modified: 21 March 2011 (DB) +// Last modified: 10 October 2011 (DB) // --------------------------------------------------------------------------- // Provides index operations for the standardized BAM index format (".bai") // *************************************************************************** -#include -#include -#include -#include -#include +#include "api/BamAlignment.h" +#include "api/internal/BamException_p.h" +#include "api/internal/BamReader_p.h" +#include "api/internal/BamStandardIndex_p.h" using namespace BamTools; using namespace BamTools::Internal; @@ -20,714 +18,724 @@ using namespace BamTools::Internal; #include #include #include -#include -#include +#include using namespace std; -BamStandardIndex::BamStandardIndex(void) - : BamIndex() - , m_dataBeginOffset(0) - , m_hasFullDataCache(false) -{ - m_isBigEndian = BamTools::SystemIsBigEndian(); -} +// ----------------------------------- +// static BamStandardIndex constants +// ----------------------------------- -BamStandardIndex::~BamStandardIndex(void) { - ClearAllData(); -} +const int BamStandardIndex::MAX_BIN = 37450; // =(8^6-1)/7+1 +const int BamStandardIndex::BAM_LIDX_SHIFT = 14; +const string BamStandardIndex::BAI_EXTENSION = ".bai"; +const char* const BamStandardIndex::BAI_MAGIC = "BAI\1"; +const int BamStandardIndex::SIZEOF_ALIGNMENTCHUNK = sizeof(uint64_t)*2; +const int BamStandardIndex::SIZEOF_BINCORE = sizeof(uint32_t) + sizeof(int32_t); +const int BamStandardIndex::SIZEOF_LINEAROFFSET = sizeof(uint64_t); -// calculate bins that overlap region -int BamStandardIndex::BinsFromRegion(const BamRegion& region, - const RefVector& references, - const bool isRightBoundSpecified, - uint16_t bins[MAX_BIN]) -{ - // get region boundaries - uint32_t begin = (unsigned int)region.LeftPosition; - uint32_t end; +// ---------------------------- +// RaiiWrapper implementation +// ---------------------------- - // if right bound specified AND left&right bounds are on same reference - // OK to use right bound position - if ( isRightBoundSpecified && ( region.LeftRefID == region.RightRefID ) ) - end = (unsigned int)region.RightPosition; +BamStandardIndex::RaiiWrapper::RaiiWrapper(void) + : IndexStream(0) + , Buffer(0) +{ } - // otherwise, use end of left bound reference as cutoff - else - end = (unsigned int)references.at(region.LeftRefID).RefLength - 1; - - // initialize list, bin '0' always a valid bin - int i = 0; - bins[i++] = 0; +BamStandardIndex::RaiiWrapper::~RaiiWrapper(void) { - // get rest of bins that contain this region - unsigned int k; - for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { bins[i++] = k; } - for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { bins[i++] = k; } - for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { bins[i++] = k; } - for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { bins[i++] = k; } - for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { bins[i++] = k; } + if ( IndexStream ) { + fclose(IndexStream); + IndexStream = 0; + } - // return number of bins stored - return i; + if ( Buffer ) { + delete[] Buffer; + Buffer = 0; + } } -// creates index data (in-memory) from @reader data -bool BamStandardIndex::Build(Internal::BamReaderPrivate* reader) { +// --------------------------------- +// BamStandardIndex implementation +// --------------------------------- - // skip if invalid reader - if ( reader == 0 ) - return false; - - // skip if reader BgzfStream is invalid or not open - BgzfStream* bgzfStream = reader->Stream(); - if ( bgzfStream == 0 || !bgzfStream->IsOpen ) - return false; - - // move reader's file pointer to beginning of alignments - reader->Rewind(); +// ctor +BamStandardIndex::BamStandardIndex(Internal::BamReaderPrivate* reader) + : BamIndex(reader) + , m_cacheMode(BamIndex::LimitedIndexCaching) + , m_bufferLength(0) +{ + m_isBigEndian = BamTools::SystemIsBigEndian(); +} - // get reference count, reserve index space - const int numReferences = reader->GetReferenceCount(); - m_indexData.clear(); - m_hasFullDataCache = false; - SetReferenceCount(numReferences); +// dtor +BamStandardIndex::~BamStandardIndex(void) { + CloseFile(); +} - // sets default constant for bin, ID, offset, coordinate variables - const uint32_t defaultValue = 0xffffffffu; +void BamStandardIndex::AdjustRegion(const BamRegion& region, uint32_t& begin, uint32_t& end) { - // bin data - uint32_t saveBin(defaultValue); - uint32_t lastBin(defaultValue); + // retrieve references from reader + const RefVector& references = m_reader->GetReferenceData(); - // reference ID data - int32_t saveRefID(defaultValue); - int32_t lastRefID(defaultValue); + // LeftPosition cannot be greater than or equal to reference length + if ( region.LeftPosition >= references.at(region.LeftRefID).RefLength ) + throw BamException("BamStandardIndex::AdjustRegion", "invalid region requested"); - // offset data - uint64_t saveOffset = bgzfStream->Tell(); - uint64_t lastOffset = saveOffset; + // set region 'begin' + begin = (unsigned int)region.LeftPosition; - // coordinate data - int32_t lastCoordinate = defaultValue; + // if right bound specified AND left&right bounds are on same reference + // OK to use right bound position as region 'end' + if ( region.isRightBoundSpecified() && ( region.LeftRefID == region.RightRefID ) ) + end = (unsigned int)region.RightPosition; - BamAlignment bAlignment; - while ( reader->LoadNextAlignment(bAlignment) ) { + // otherwise, set region 'end' to last reference base + else end = (unsigned int)references.at(region.LeftRefID).RefLength; +} - // change of chromosome, save ID, reset bin - if ( lastRefID != bAlignment.RefID ) { - lastRefID = bAlignment.RefID; - lastBin = defaultValue; - } +// [begin, end) +void BamStandardIndex::CalculateCandidateBins(const uint32_t& begin, + const uint32_t& end, + set& candidateBins) +{ + // initialize list, bin '0' is always a valid bin + candidateBins.insert(0); - // if lastCoordinate greater than BAM position - file not sorted properly - else if ( lastCoordinate > bAlignment.Position ) { - fprintf(stderr, "BamStandardIndex ERROR: file not properly sorted:\n"); - fprintf(stderr, "Alignment %s : %d > %d on reference (id = %d)", - bAlignment.Name.c_str(), lastCoordinate, bAlignment.Position, bAlignment.RefID); - exit(1); - } + // get rest of bins that contain this region + unsigned int k; + for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { candidateBins.insert(k); } + for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { candidateBins.insert(k); } + for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { candidateBins.insert(k); } + for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { candidateBins.insert(k); } + for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { candidateBins.insert(k); } +} - // if valid reference && BAM bin spans some minimum cutoff (smaller bin ids span larger regions) - if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) { +void BamStandardIndex::CalculateCandidateOffsets(const BaiReferenceSummary& refSummary, + const uint64_t& minOffset, + set& candidateBins, + vector& offsets) +{ + // seek to first bin + Seek(refSummary.FirstBinFilePosition, SEEK_SET); - // save linear offset entry (matched to BAM entry refID) - BamStandardIndexData::iterator indexIter = m_indexData.find(bAlignment.RefID); - if ( indexIter == m_indexData.end() ) return false; // error - ReferenceIndex& refIndex = (*indexIter).second; - LinearOffsetVector& offsets = refIndex.Offsets; - SaveLinearOffset(offsets, bAlignment, lastOffset); - } + // iterate over reference bins + uint32_t binId; + int32_t numAlignmentChunks; + set::iterator candidateBinIter; + for ( int i = 0; i < refSummary.NumBins; ++i ) { - // if current BamAlignment bin != lastBin, "then possibly write the binning index" - if ( bAlignment.Bin != lastBin ) { + // read bin contents (if successful, alignment chunks are now in m_buffer) + ReadBinIntoBuffer(binId, numAlignmentChunks); - // if not first time through - if ( saveBin != defaultValue ) { + // see if bin is a 'candidate bin' + candidateBinIter = candidateBins.find(binId); - // save Bam bin entry - BamStandardIndexData::iterator indexIter = m_indexData.find(saveRefID); - if ( indexIter == m_indexData.end() ) return false; // error - ReferenceIndex& refIndex = (*indexIter).second; - BamBinMap& binMap = refIndex.Bins; - SaveBinEntry(binMap, saveBin, saveOffset, lastOffset); - } + // if not, move on to next bin + if ( candidateBinIter == candidateBins.end() ) + continue; - // update saveOffset - saveOffset = lastOffset; + // otherwise, check bin's contents against for overlap + else { - // update bin values - saveBin = bAlignment.Bin; - lastBin = bAlignment.Bin; + size_t offset = 0; + uint64_t chunkStart; + uint64_t chunkStop; - // update saveRefID - saveRefID = bAlignment.RefID; + // iterate over alignment chunks + for ( int j = 0; j < numAlignmentChunks; ++j ) { - // if invalid RefID, break out - if ( saveRefID < 0 ) break; - } + // read chunk start & stop from buffer + memcpy((char*)&chunkStart, Resources.Buffer+offset, sizeof(uint64_t)); + offset += sizeof(uint64_t); + memcpy((char*)&chunkStop, Resources.Buffer+offset, sizeof(uint64_t)); + offset += sizeof(uint64_t); - // make sure that current file pointer is beyond lastOffset - if ( bgzfStream->Tell() <= (int64_t)lastOffset ) { - fprintf(stderr, "BamStandardIndex ERROR: could not build index - calculating offsets failed.\n"); - exit(1); - } + // swap endian-ness if necessary + if ( m_isBigEndian ) { + SwapEndian_64(chunkStart); + SwapEndian_64(chunkStop); + } - // update lastOffset - lastOffset = bgzfStream->Tell(); + // store alignment chunk's start offset + // if its stop offset is larger than our 'minOffset' + if ( chunkStop >= minOffset ) + offsets.push_back(chunkStart); + } - // update lastCoordinate - lastCoordinate = bAlignment.Position; - } + // 'pop' bin ID from candidate bins set + candidateBins.erase(candidateBinIter); - // save any leftover BAM data (as long as refID is valid) - if ( saveRefID >= 0 ) { - // save Bam bin entry - BamStandardIndexData::iterator indexIter = m_indexData.find(saveRefID); - if ( indexIter == m_indexData.end() ) return false; // error - ReferenceIndex& refIndex = (*indexIter).second; - BamBinMap& binMap = refIndex.Bins; - SaveBinEntry(binMap, saveBin, saveOffset, lastOffset); + // quit if no more candidates + if ( candidateBins.empty() ) + break; + } } +} - // simplify index by merging chunks - MergeChunks(); - - // iterate through references in index - // sort offsets in linear offset vector - BamStandardIndexData::iterator indexIter = m_indexData.begin(); - BamStandardIndexData::iterator indexEnd = m_indexData.end(); - for ( int i = 0; indexIter != indexEnd; ++indexIter, ++i ) { - - // get reference index data - ReferenceIndex& refIndex = (*indexIter).second; - LinearOffsetVector& offsets = refIndex.Offsets; +uint64_t BamStandardIndex::CalculateMinOffset(const BaiReferenceSummary& refSummary, + const uint32_t& begin) +{ + // if no linear offsets exist, return 0 + if ( refSummary.NumLinearOffsets == 0 ) + return 0; + + // if 'begin' starts beyond last linear offset, use the last linear offset as minimum + // else use the offset corresponding to the requested start position + const int shiftedBegin = begin>>BamStandardIndex::BAM_LIDX_SHIFT; + if ( shiftedBegin >= refSummary.NumLinearOffsets ) + return LookupLinearOffset( refSummary, refSummary.NumLinearOffsets-1 ); + else + return LookupLinearOffset( refSummary, shiftedBegin ); +} - // sort linear offsets - sort(offsets.begin(), offsets.end()); +void BamStandardIndex::CheckBufferSize(char*& buffer, + unsigned int& bufferLength, + const unsigned int& requestedBytes) +{ + try { + if ( requestedBytes > bufferLength ) { + bufferLength = requestedBytes + 10; + delete[] buffer; + buffer = new char[bufferLength]; + } + } catch ( std::bad_alloc& ) { + stringstream s(""); + s << "out of memory when allocating " << requestedBytes << " bytes"; + throw BamException("BamStandardIndex::CheckBufferSize", s.str()); } +} - // rewind reader's file pointer to beginning of alignments, return success/fail - return reader->Rewind(); +void BamStandardIndex::CheckBufferSize(unsigned char*& buffer, + unsigned int& bufferLength, + const unsigned int& requestedBytes) +{ + try { + if ( requestedBytes > bufferLength ) { + bufferLength = requestedBytes + 10; + delete[] buffer; + buffer = new unsigned char[bufferLength]; + } + } catch ( std::bad_alloc& ) { + stringstream s(""); + s << "out of memory when allocating " << requestedBytes << " bytes"; + throw BamException("BamStandardIndex::CheckBufferSize", s.str()); + } } -// check index file magic number, return true if OK -bool BamStandardIndex::CheckMagicNumber(void) { +void BamStandardIndex::CheckMagicNumber(void) { - // read in magic number + // check 'magic number' to see if file is BAI index char magic[4]; - size_t elementsRead = fread(magic, sizeof(char), 4, m_indexStream); + const size_t elementsRead = fread(magic, sizeof(char), 4, Resources.IndexStream); + if ( elementsRead != 4 ) + throw BamException("BamStandardIndex::CheckMagicNumber", "could not read BAI magic number"); // compare to expected value - if ( strncmp(magic, "BAI\1", 4) != 0 ) { - fprintf(stderr, "BamStandardIndex ERROR: could not load index file - invalid magic number.\n"); - fclose(m_indexStream); - return false; - } - - // return success/failure of load - return (elementsRead == 4); + if ( strncmp(magic, BamStandardIndex::BAI_MAGIC, 4) != 0 ) + throw BamException("BamStandardIndex::CheckMagicNumber", "invalid BAI magic number"); } -// clear all current index offset data in memory -void BamStandardIndex::ClearAllData(void) { - BamStandardIndexData::const_iterator indexIter = m_indexData.begin(); - BamStandardIndexData::const_iterator indexEnd = m_indexData.end(); - for ( ; indexIter != indexEnd; ++indexIter ) { - const int& refId = (*indexIter).first; - ClearReferenceOffsets(refId); - } +void BamStandardIndex::ClearReferenceEntry(BaiReferenceEntry& refEntry) { + refEntry.ID = -1; + refEntry.Bins.clear(); + refEntry.LinearOffsets.clear(); } -// clear all index offset data for desired reference -void BamStandardIndex::ClearReferenceOffsets(const int& refId) { +void BamStandardIndex::CloseFile(void) { - // look up refId, skip if not found - BamStandardIndexData::iterator indexIter = m_indexData.find(refId); - if ( indexIter == m_indexData.end() ) return ; + // close file stream + if ( IsFileOpen() ) { + fclose(Resources.IndexStream); + Resources.IndexStream = 0; + } - // clear reference data - ReferenceIndex& refEntry = (*indexIter).second; - refEntry.Bins.clear(); - refEntry.Offsets.clear(); + // clear index file summary data + m_indexFileSummary.clear(); - // set flag - m_hasFullDataCache = false; + // clean up I/O buffer + delete[] Resources.Buffer; + Resources.Buffer = 0; + m_bufferLength = 0; } -// return file position after header metadata -off_t BamStandardIndex::DataBeginOffset(void) const { - return m_dataBeginOffset; -} +// builds index from associated BAM file & writes out to index file +bool BamStandardIndex::Create(void) { -// calculates offset(s) for a given region -bool BamStandardIndex::GetOffsets(const BamRegion& region, - const RefVector& references, - const bool isRightBoundSpecified, - vector& offsets, - bool* hasAlignmentsInRegion) -{ - // return false if leftBound refID is not found in index data - if ( m_indexData.find(region.LeftRefID) == m_indexData.end() ) + // skip if BamReader is invalid or not open + if ( m_reader == 0 || !m_reader->IsOpen() ) { + SetErrorString("BamStandardIndex::Create", "could not create index: reader is not open"); return false; - - // load index data for region if not already cached - if ( !IsDataLoaded(region.LeftRefID) ) { - bool loadedOk = true; - loadedOk &= SkipToReference(region.LeftRefID); - loadedOk &= LoadReference(region.LeftRefID); - if ( !loadedOk ) return false; } - // calculate which bins overlap this region - uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2); - int numBins = BinsFromRegion(region, references, isRightBoundSpecified, bins); - - // get bins for this reference - BamStandardIndexData::const_iterator indexIter = m_indexData.find(region.LeftRefID); - if ( indexIter == m_indexData.end() ) return false; // error - const ReferenceIndex& refIndex = (*indexIter).second; - const BamBinMap& binMap = refIndex.Bins; - - // get minimum offset to consider - const LinearOffsetVector& linearOffsets = refIndex.Offsets; - const uint64_t minOffset = ( (unsigned int)(region.LeftPosition>>BAM_LIDX_SHIFT) >= linearOffsets.size() ) - ? 0 : linearOffsets.at(region.LeftPosition>>BAM_LIDX_SHIFT); - - // store all alignment 'chunk' starts (file offsets) for bins in this region - for ( int i = 0; i < numBins; ++i ) { - - const uint16_t binKey = bins[i]; - map::const_iterator binIter = binMap.find(binKey); - if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) { - - // iterate over chunks - const ChunkVector& chunks = (*binIter).second; - std::vector::const_iterator chunksIter = chunks.begin(); - std::vector::const_iterator chunksEnd = chunks.end(); - for ( ; chunksIter != chunksEnd; ++chunksIter) { - - // if valid chunk found, store its file offset - const Chunk& chunk = (*chunksIter); - if ( chunk.Stop > minOffset ) - offsets.push_back( chunk.Start ); - } - } + // rewind BamReader + if ( !m_reader->Rewind() ) { + const string readerError = m_reader->GetErrorString(); + const string message = "could not create index: \n\t" + readerError; + SetErrorString("BamStandardIndex::Create", message); + return false; } - // clean up memory - free(bins); - - // sort the offsets before returning - sort(offsets.begin(), offsets.end()); + try { + + // open new index file (read & write) + string indexFilename = m_reader->Filename() + Extension(); + OpenFile(indexFilename, "w+b"); + + // initialize BaiFileSummary with number of references + const int& numReferences = m_reader->GetReferenceCount(); + ReserveForSummary(numReferences); + + // initialize output file + WriteHeader(); + + // set up bin, ID, offset, & coordinate markers + const uint32_t defaultValue = 0xffffffffu; + uint32_t currentBin = defaultValue; + uint32_t lastBin = defaultValue; + int32_t currentRefID = defaultValue; + int32_t lastRefID = defaultValue; + uint64_t currentOffset = (uint64_t)m_reader->Tell(); + uint64_t lastOffset = currentOffset; + int32_t lastPosition = defaultValue; + + // iterate through alignments in BAM file + BamAlignment al; + BaiReferenceEntry refEntry; + while ( m_reader->LoadNextAlignment(al) ) { + + // changed to new reference + if ( lastRefID != al.RefID ) { + + // if not first reference, save previous reference data + if ( lastRefID != (int32_t)defaultValue ) { + + SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset); + WriteReferenceEntry(refEntry); + ClearReferenceEntry(refEntry); + + // write any empty references between (but *NOT* including) lastRefID & al.RefID + for ( int i = lastRefID+1; i < al.RefID; ++i ) { + BaiReferenceEntry emptyEntry(i); + WriteReferenceEntry(emptyEntry); + } + + // update bin markers + currentOffset = lastOffset; + currentBin = al.Bin; + lastBin = al.Bin; + currentRefID = al.RefID; + } - // set flag & return success - *hasAlignmentsInRegion = (offsets.size() != 0 ); + // otherwise, this is first pass + // be sure to write any empty references up to (but *NOT* including) current RefID + else { + for ( int i = 0; i < al.RefID; ++i ) { + BaiReferenceEntry emptyEntry(i); + WriteReferenceEntry(emptyEntry); + } + } - // if cache mode set to none, dump the data we just loaded - if ( m_cacheMode == BamIndex::NoIndexCaching ) - ClearReferenceOffsets(region.LeftRefID); + // update reference markers + refEntry.ID = al.RefID; + lastRefID = al.RefID; + lastBin = defaultValue; + } - // return succes - return true; -} + // if lastPosition greater than current alignment position - file not sorted properly + else if ( lastPosition > al.Position ) { + stringstream s(""); + s << "BAM file is not properly sorted by coordinate" << endl + << "Current alignment position: " << al.Position + << " < previous alignment position: " << lastPosition + << " on reference ID: " << al.RefID << endl; + SetErrorString("BamStandardIndex::Create", s.str()); + return false; + } -// returns whether reference has alignments or no -bool BamStandardIndex::HasAlignments(const int& refId) const { - BamStandardIndexData::const_iterator indexIter = m_indexData.find(refId); - if ( indexIter == m_indexData.end() ) return false; // error - const ReferenceIndex& refEntry = (*indexIter).second; - return refEntry.HasAlignments; -} + // if alignment's ref ID is valid & its bin is not a 'leaf' + if ( (al.RefID >= 0) && (al.Bin < 4681) ) + SaveLinearOffsetEntry(refEntry.LinearOffsets, al.Position, al.GetEndPosition(), lastOffset); -// return true if all index data is cached -bool BamStandardIndex::HasFullDataCache(void) const { - return m_hasFullDataCache; -} + // changed to new BAI bin + if ( al.Bin != lastBin ) { -// returns true if index cache has data for desired reference -bool BamStandardIndex::IsDataLoaded(const int& refId) const { + // if not first bin on reference, save previous bin data + if ( currentBin != defaultValue ) + SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset); - // look up refId, return false if not found - BamStandardIndexData::const_iterator indexIter = m_indexData.find(refId); - if ( indexIter == m_indexData.end() ) return false; + // update markers + currentOffset = lastOffset; + currentBin = al.Bin; + lastBin = al.Bin; + currentRefID = al.RefID; - // see if reference has alignments - // if not, it's not a problem to have no offset data - const ReferenceIndex& refEntry = (*indexIter).second; - if ( !refEntry.HasAlignments ) return true; + // if invalid RefID, break out + if ( currentRefID < 0 ) + break; + } - // return whether bin map contains data - return ( !refEntry.Bins.empty() ); -} + // make sure that current file pointer is beyond lastOffset + if ( m_reader->Tell() <= (int64_t)lastOffset ) { + SetErrorString("BamStandardIndex::Create", "calculating offsets failed"); + return false; + } -// attempts to use index to jump to region; returns success/fail -bool BamStandardIndex::Jump(Internal::BamReaderPrivate* reader, - const BamTools::BamRegion& region, - bool *hasAlignmentsInRegion) -{ - // skip if invalid reader - if ( reader == 0 ) - return false; + // update lastOffset & lastPosition + lastOffset = m_reader->Tell(); + lastPosition = al.Position; + } - // skip if reader BgzfStream is invalid or not open - BgzfStream* bgzfStream = reader->Stream(); - if ( bgzfStream == 0 || !bgzfStream->IsOpen ) - return false; + // after finishing alignments, if any data was read, check: + if ( currentRefID >= 0 ) { - // retrieve references from reader - const RefVector references = reader->GetReferenceData(); + // store last alignment chunk to its bin, then write last reference entry with data + SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset); + WriteReferenceEntry(refEntry); - // make sure left-bound position is valid - if ( region.LeftPosition > references.at(region.LeftRefID).RefLength ) - return false; + // then write any empty references remaining at end of file + for ( int i = currentRefID+1; i < numReferences; ++i ) { + BaiReferenceEntry emptyEntry(i); + WriteReferenceEntry(emptyEntry); + } + } - // calculate offsets for this region - // if failed, print message, set flag, and return failure - vector offsets; - if ( !GetOffsets(region, references, region.isRightBoundSpecified(), offsets, hasAlignmentsInRegion) ) { - fprintf(stderr, "BamStandardIndex ERROR: could not jump - unable to calculate offset candidates for specified region.\n"); - *hasAlignmentsInRegion = false; + } catch ( BamException& e) { + m_errorString = e.what(); return false; } - // iterate through offsets - BamAlignment alignment; - bool result = true; - for ( vector::const_iterator o = offsets.begin(); o != offsets.end(); ++o) { - - // attempt seek & load first available alignment - // set flag to true if data exists - result &= bgzfStream->Seek(*o); - *hasAlignmentsInRegion = reader->GetNextAlignmentCore(alignment); - - // if this alignment corresponds to desired position - // return success of seeking back to the offset before the 'current offset' (to cover overlaps) - if ( ((alignment.RefID == region.LeftRefID) && - ((alignment.Position + alignment.Length) > region.LeftPosition)) || - (alignment.RefID > region.LeftRefID) ) - { - if ( o != offsets.begin() ) --o; - return bgzfStream->Seek(*o); - } - } - - // if error in jumping, print message & set flag - if ( !result ) { - fprintf(stderr, "BamStandardIndex ERROR: could not jump - unable to determine correct offset for specified region.\n"); - *hasAlignmentsInRegion = false; + // rewind BamReader + if ( !m_reader->Rewind() ) { + const string readerError = m_reader->GetErrorString(); + const string message = "could not create index: \n\t" + readerError; + SetErrorString("BamStandardIndex::Create", message); + return false; } - // return success/failure - return result; -} - -// clears index data from all references except the first -void BamStandardIndex::KeepOnlyFirstReferenceOffsets(void) { - BamStandardIndexData::const_iterator indexBegin = m_indexData.begin(); - KeepOnlyReferenceOffsets((*indexBegin).first); + // return success + return true; } -// clears index data from all references except the one specified -void BamStandardIndex::KeepOnlyReferenceOffsets(const int& refId) { - BamStandardIndexData::iterator mapIter = m_indexData.begin(); - BamStandardIndexData::iterator mapEnd = m_indexData.end(); - for ( ; mapIter != mapEnd; ++mapIter ) { - const int entryRefId = (*mapIter).first; - if ( entryRefId != refId ) - ClearReferenceOffsets(entryRefId); - } +// returns format's file extension +const string BamStandardIndex::Extension(void) { + return BamStandardIndex::BAI_EXTENSION; } -bool BamStandardIndex::LoadAllReferences(bool saveData) { +void BamStandardIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) { - // skip if data already loaded - if ( m_hasFullDataCache ) return true; + // cannot calculate offsets if unknown/invalid reference ID requested + if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() ) + throw BamException("BamStandardIndex::GetOffset", "invalid reference ID requested"); - // get number of reference sequences - uint32_t numReferences; - if ( !LoadReferenceCount((int&)numReferences) ) - return false; + // retrieve index summary for left bound reference + const BaiReferenceSummary& refSummary = m_indexFileSummary.at(region.LeftRefID); - // iterate over reference entries - bool loadedOk = true; - for ( int i = 0; i < (int)numReferences; ++i ) - loadedOk &= LoadReference(i, saveData); + // set up region boundaries based on actual BamReader data + uint32_t begin; + uint32_t end; + AdjustRegion(region, begin, end); - // set flag - if ( loadedOk && saveData ) - m_hasFullDataCache = true; + // retrieve all candidate bin IDs for region + set candidateBins; + CalculateCandidateBins(begin, end, candidateBins); - // return success/failure of loading references - return loadedOk; -} + // use reference's linear offsets to calculate the minimum offset + // that must be considered to find overlap + const uint64_t& minOffset = CalculateMinOffset(refSummary, begin); -// load header data from index file, return true if loaded OK -bool BamStandardIndex::LoadHeader(void) { + // attempt to use reference summary, minOffset, & candidateBins to calculate offsets + // no data should not be error, just bail + vector offsets; + CalculateCandidateOffsets(refSummary, minOffset, candidateBins, offsets); + if ( offsets.empty() ) + return; + + // ensure that offsets are sorted before processing + sort( offsets.begin(), offsets.end() ); + + // binary search for an overlapping block (may not be first one though) + BamAlignment al; + typedef vector::const_iterator OffsetConstIterator; + OffsetConstIterator offsetFirst = offsets.begin(); + OffsetConstIterator offsetIter = offsetFirst; + OffsetConstIterator offsetLast = offsets.end(); + iterator_traits::difference_type count = distance(offsetFirst, offsetLast); + iterator_traits::difference_type step; + while ( count > 0 ) { + offsetIter = offsetFirst; + step = count/2; + advance(offsetIter, step); + + // attempt seek to candidate offset + const int64_t& candidateOffset = (*offsetIter); + if ( !m_reader->Seek(candidateOffset) ) { + const string readerError = m_reader->GetErrorString(); + const string message = "could not seek in BAM file: \n\t" + readerError; + throw BamException("BamToolsIndex::GetOffset", message); + } - bool loadedOk = CheckMagicNumber(); + // load first available alignment, setting flag to true if data exists + *hasAlignmentsInRegion = m_reader->LoadNextAlignment(al); - // store offset of beginning of data - m_dataBeginOffset = ftell64(m_indexStream); + // check alignment against region + if ( al.GetEndPosition() <= region.LeftPosition ) { + offsetFirst = ++offsetIter; + count -= step+1; + } else count = step; + } - // return success/failure of load - return loadedOk; + // step back to the offset before the 'current offset' (to make sure we cover overlaps) + if ( offsetIter != offsets.begin() ) + --offsetIter; + offset = (*offsetIter); } -// load a single index bin entry from file, return true if loaded OK -// @saveData - save data in memory if true, just read & discard if false -bool BamStandardIndex::LoadBin(ReferenceIndex& refEntry, bool saveData) { - - size_t elementsRead = 0; - - // get bin ID - uint32_t binId; - elementsRead += fread(&binId, sizeof(binId), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_32(binId); - - // load alignment chunks for this bin - ChunkVector chunks; - bool chunksOk = LoadChunks(chunks, saveData); - - // store bin entry - if ( chunksOk && saveData ) - refEntry.Bins.insert(pair(binId, chunks)); - - // return success/failure of load - return ( (elementsRead == 1) && chunksOk ); +// returns whether reference has alignments or no +bool BamStandardIndex::HasAlignments(const int& referenceID) const { + if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() ) + return false; + const BaiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID); + return ( refSummary.NumBins > 0 ); } -bool BamStandardIndex::LoadBins(ReferenceIndex& refEntry, bool saveData) { - - size_t elementsRead = 0; - - // get number of bins - int32_t numBins; - elementsRead += fread(&numBins, sizeof(numBins), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_32(numBins); - - // set flag - refEntry.HasAlignments = ( numBins != 0 ); - - // iterate over bins - bool binsOk = true; - for ( int i = 0; i < numBins; ++i ) - binsOk &= LoadBin(refEntry, saveData); - - // return success/failure of load - return ( (elementsRead == 1) && binsOk ); +bool BamStandardIndex::IsFileOpen(void) const { + return ( Resources.IndexStream != 0 ); } -// load a single index bin entry from file, return true if loaded OK -// @saveData - save data in memory if true, just read & discard if false -bool BamStandardIndex::LoadChunk(ChunkVector& chunks, bool saveData) { +// attempts to use index data to jump to @region, returns success/fail +// a "successful" jump indicates no error, but not whether this region has data +// * thus, the method sets a flag to indicate whether there are alignments +// available after the jump position +bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) { - size_t elementsRead = 0; + // clear out flag + *hasAlignmentsInRegion = false; - // read in chunk data - uint64_t start; - uint64_t stop; - elementsRead += fread(&start, sizeof(start), 1, m_indexStream); - elementsRead += fread(&stop, sizeof(stop), 1, m_indexStream); + // skip if invalid reader or not open + if ( m_reader == 0 || !m_reader->IsOpen() ) { + SetErrorString("BamStandardIndex::Jump", "could not jump: reader is not open"); + return false; + } - // swap endian-ness if necessary - if ( m_isBigEndian ) { - SwapEndian_64(start); - SwapEndian_64(stop); + // calculate nearest offset to jump to + int64_t offset; + try { + GetOffset(region, offset, hasAlignmentsInRegion); + } catch ( BamException& e ) { + m_errorString = e.what(); + return false; } - // save data if requested - if ( saveData ) chunks.push_back( Chunk(start, stop) ); + // if region has alignments, return success/fail of seeking there + if ( *hasAlignmentsInRegion ) + return m_reader->Seek(offset); - // return success/failure of load - return ( elementsRead == 2 ); + // otherwise, simply return true (but hasAlignmentsInRegion flag has been set to false) + // (this is OK, BamReader will check this flag before trying to load data) + return true; } -bool BamStandardIndex::LoadChunks(ChunkVector& chunks, bool saveData) { +// loads existing data from file into memory +bool BamStandardIndex::Load(const std::string& filename) { - size_t elementsRead = 0; + try { - // read in number of chunks - uint32_t numChunks; - elementsRead += fread(&numChunks, sizeof(numChunks), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_32(numChunks); + // attempt to open file (read-only) + OpenFile(filename, "rb"); - // initialize space for chunks if we're storing this data - if ( saveData ) chunks.reserve(numChunks); + // validate format + CheckMagicNumber(); - // iterate over chunks - bool chunksOk = true; - for ( int i = 0; i < (int)numChunks; ++i ) - chunksOk &= LoadChunk(chunks, saveData); + // load in-memory summary of index data + SummarizeIndexFile(); - // sort chunk vector - sort( chunks.begin(), chunks.end(), ChunkLessThan ); + // return success + return true; - // return success/failure of load - return ( (elementsRead == 1) && chunksOk ); + } catch ( BamException& e ) { + m_errorString = e.what(); + return false; + } } -// load a single index linear offset entry from file, return true if loaded OK -// @saveData - save data in memory if true, just read & discard if false -bool BamStandardIndex::LoadLinearOffsets(ReferenceIndex& refEntry, bool saveData) { +uint64_t BamStandardIndex::LookupLinearOffset(const BaiReferenceSummary& refSummary, const int& index) { - size_t elementsRead = 0; + // attempt seek to proper index file position + const int64_t linearOffsetFilePosition = (int64_t)refSummary.FirstLinearOffsetFilePosition + + index*BamStandardIndex::SIZEOF_LINEAROFFSET; + Seek(linearOffsetFilePosition, SEEK_SET); - // read in number of linear offsets - int32_t numLinearOffsets; - elementsRead += fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets); + // read linear offset from BAI file + uint64_t linearOffset; + ReadLinearOffset(linearOffset); + return linearOffset; +} - // set up destination vector (if we're saving the data) - LinearOffsetVector linearOffsets; - if ( saveData ) linearOffsets.reserve(numLinearOffsets); +void BamStandardIndex::MergeAlignmentChunks(BaiAlignmentChunkVector& chunks) { - // iterate over linear offsets - uint64_t linearOffset; - for ( int i = 0; i < numLinearOffsets; ++i ) { - elementsRead += fread(&linearOffset, sizeof(linearOffset), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_64(linearOffset); - if ( saveData ) linearOffsets.push_back(linearOffset); - } + // skip if chunks are empty, nothing to merge + if ( chunks.empty() ) + return; - // sort linear offsets - sort ( linearOffsets.begin(), linearOffsets.end() ); + // set up merged alignment chunk container + BaiAlignmentChunkVector mergedChunks; + mergedChunks.push_back( chunks[0] ); - // save in reference index entry if desired - if ( saveData ) refEntry.Offsets = linearOffsets; + // iterate over chunks + int i = 0; + BaiAlignmentChunkVector::iterator chunkIter = chunks.begin(); + BaiAlignmentChunkVector::iterator chunkEnd = chunks.end(); + for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) { - // return success/failure of load - return ( elementsRead == (size_t)(numLinearOffsets + 1) ); -} + // get 'currentMergeChunk' based on numeric index + BaiAlignmentChunk& currentMergeChunk = mergedChunks[i]; -bool BamStandardIndex::LoadFirstReference(bool saveData) { - BamStandardIndexData::const_iterator indexBegin = m_indexData.begin(); - return LoadReference((*indexBegin).first, saveData); -} + // get sourceChunk based on source vector iterator + BaiAlignmentChunk& sourceChunk = (*chunkIter); -// load a single reference from file, return true if loaded OK -// @saveData - save data in memory if true, just read & discard if false -bool BamStandardIndex::LoadReference(const int& refId, bool saveData) { + // if currentMergeChunk ends where sourceChunk starts, then merge the two + if ( currentMergeChunk.Stop>>16 == sourceChunk.Start>>16 ) + currentMergeChunk.Stop = sourceChunk.Stop; - // look up refId - BamStandardIndexData::iterator indexIter = m_indexData.find(refId); + // otherwise + else { + // append sourceChunk after currentMergeChunk + mergedChunks.push_back(sourceChunk); - // if reference not previously loaded, create new entry - if ( indexIter == m_indexData.end() ) { - ReferenceIndex newEntry; - newEntry.HasAlignments = false; - m_indexData.insert( pair(refId, newEntry) ); + // update i, so the next iteration will consider the + // recently-appended sourceChunk as new mergeChunk candidate + ++i; + } } - // load reference data - indexIter = m_indexData.find(refId); - ReferenceIndex& entry = (*indexIter).second; - bool loadedOk = true; - loadedOk &= LoadBins(entry, saveData); - loadedOk &= LoadLinearOffsets(entry, saveData); - return loadedOk; + // saved newly-merged chunks into (parameter) chunks + chunks = mergedChunks; } -// loads number of references, return true if loaded OK -bool BamStandardIndex::LoadReferenceCount(int& numReferences) { +void BamStandardIndex::OpenFile(const std::string& filename, const char* mode) { - size_t elementsRead = 0; + // make sure any previous index file is closed + CloseFile(); - // read reference count - elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream); - if ( m_isBigEndian ) SwapEndian_32(numReferences); + // attempt to open file + Resources.IndexStream = fopen(filename.c_str(), mode); + if ( !IsFileOpen() ) { + const string message = string("could not open file: ") + filename; + throw BamException("BamStandardIndex::OpenFile", message); + } +} - // return success/failure of load - return ( elementsRead == 1 ); +void BamStandardIndex::ReadBinID(uint32_t& binId) { + const size_t elementsRead = fread(&binId, sizeof(binId), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_32(binId); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadBinID", "could not read BAI bin ID"); } -// merges 'alignment chunks' in BAM bin (used for index building) -void BamStandardIndex::MergeChunks(void) { +void BamStandardIndex::ReadBinIntoBuffer(uint32_t& binId, int32_t& numAlignmentChunks) { - // iterate over reference enties - BamStandardIndexData::iterator indexIter = m_indexData.begin(); - BamStandardIndexData::iterator indexEnd = m_indexData.end(); - for ( ; indexIter != indexEnd; ++indexIter ) { + // read bin header + ReadBinID(binId); + ReadNumAlignmentChunks(numAlignmentChunks); - // get BAM bin map for this reference - ReferenceIndex& refIndex = (*indexIter).second; - BamBinMap& bamBinMap = refIndex.Bins; + // read bin contents + const unsigned int bytesRequested = numAlignmentChunks*BamStandardIndex::SIZEOF_ALIGNMENTCHUNK; + ReadIntoBuffer(bytesRequested); +} - // iterate over BAM bins - BamBinMap::iterator binIter = bamBinMap.begin(); - BamBinMap::iterator binEnd = bamBinMap.end(); - for ( ; binIter != binEnd; ++binIter ) { +void BamStandardIndex::ReadIntoBuffer(const unsigned int& bytesRequested) { - // get chunk vector for this bin - ChunkVector& binChunks = (*binIter).second; - if ( binChunks.size() == 0 ) continue; + // ensure that our buffer is big enough for request + BamStandardIndex::CheckBufferSize(Resources.Buffer, m_bufferLength, bytesRequested); - ChunkVector mergedChunks; - mergedChunks.push_back( binChunks[0] ); + // read from BAI file stream + const size_t bytesRead = fread( Resources.Buffer, sizeof(char), bytesRequested, Resources.IndexStream ); + if ( bytesRead != (size_t)bytesRequested ) { + stringstream s(""); + s << "expected to read: " << bytesRequested << " bytes, " + << "but instead read: " << bytesRead; + throw BamException("BamStandardIndex::ReadIntoBuffer", s.str()); + } +} - // iterate over chunks - int i = 0; - ChunkVector::iterator chunkIter = binChunks.begin(); - ChunkVector::iterator chunkEnd = binChunks.end(); - for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) { +void BamStandardIndex::ReadLinearOffset(uint64_t& linearOffset) { + const size_t elementsRead = fread(&linearOffset, sizeof(linearOffset), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_64(linearOffset); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadLinearOffset", "could not read BAI linear offset"); +} - // get 'currentChunk' based on numeric index - Chunk& currentChunk = mergedChunks[i]; +void BamStandardIndex::ReadNumAlignmentChunks(int& numAlignmentChunks) { + const size_t elementsRead = fread(&numAlignmentChunks, sizeof(numAlignmentChunks), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_32(numAlignmentChunks); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadNumAlignmentChunks", "could not read BAI chunk count"); +} - // get iteratorChunk based on vector iterator - Chunk& iteratorChunk = (*chunkIter); +void BamStandardIndex::ReadNumBins(int& numBins) { + const size_t elementsRead = fread(&numBins, sizeof(numBins), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_32(numBins); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadNumBins", "could not read BAI bin count"); +} - // if chunk ends where (iterator) chunk starts, then merge - if ( currentChunk.Stop>>16 == iteratorChunk.Start>>16 ) - currentChunk.Stop = iteratorChunk.Stop; +void BamStandardIndex::ReadNumLinearOffsets(int& numLinearOffsets) { + const size_t elementsRead = fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadNumAlignmentChunks", "could not read BAI linear offset count"); +} - // otherwise - else { - // set currentChunk + 1 to iteratorChunk - mergedChunks.push_back(iteratorChunk); - ++i; - } - } +void BamStandardIndex::ReadNumReferences(int& numReferences) { + const size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, Resources.IndexStream); + if ( m_isBigEndian ) SwapEndian_32(numReferences); + if ( elementsRead != 1 ) + throw BamException("BamStandardIndex::ReadNumReferences", "could not read reference count"); +} - // saved merged chunk vector - (*binIter).second = mergedChunks; - } - } +void BamStandardIndex::ReserveForSummary(const int& numReferences) { + m_indexFileSummary.clear(); + m_indexFileSummary.assign( numReferences, BaiReferenceSummary() ); } -// saves BAM bin entry for index -void BamStandardIndex::SaveBinEntry(BamBinMap& binMap, - const uint32_t& saveBin, - const uint64_t& saveOffset, - const uint64_t& lastOffset) +void BamStandardIndex::SaveAlignmentChunkToBin(BaiBinMap& binMap, + const uint32_t& currentBin, + const uint64_t& currentOffset, + const uint64_t& lastOffset) { - // look up saveBin - BamBinMap::iterator binIter = binMap.find(saveBin); - - // create new chunk - Chunk newChunk(saveOffset, lastOffset); + // create new alignment chunk + BaiAlignmentChunk newChunk(currentOffset, lastOffset); - // if entry doesn't exist + // if no entry exists yet for this bin, create one and store alignment chunk + BaiBinMap::iterator binIter = binMap.find(currentBin); if ( binIter == binMap.end() ) { - ChunkVector newChunks; + BaiAlignmentChunkVector newChunks; newChunks.push_back(newChunk); - binMap.insert( pair(saveBin, newChunks)); + binMap.insert( pair(currentBin, newChunks)); } - // otherwise + // otherwise, just append alignment chunk else { - ChunkVector& binChunks = (*binIter).second; + BaiAlignmentChunkVector& binChunks = (*binIter).second; binChunks.push_back( newChunk ); } } -// saves linear offset entry for index -void BamStandardIndex::SaveLinearOffset(LinearOffsetVector& offsets, - const BamAlignment& bAlignment, - const uint64_t& lastOffset) +void BamStandardIndex::SaveBinsSummary(const int& refId, const int& numBins) { + BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId); + refSummary.NumBins = numBins; + refSummary.FirstBinFilePosition = Tell(); +} + +void BamStandardIndex::SaveLinearOffsetEntry(BaiLinearOffsetVector& offsets, + const int& alignmentStartPosition, + const int& alignmentStopPosition, + const uint64_t& lastOffset) { // get converted offsets - int beginOffset = bAlignment.Position >> BAM_LIDX_SHIFT; - int endOffset = (bAlignment.GetEndPosition() - 1) >> BAM_LIDX_SHIFT; + const int beginOffset = alignmentStartPosition >> BamStandardIndex::BAM_LIDX_SHIFT; + const int endOffset = (alignmentStopPosition - 1) >> BamStandardIndex::BAM_LIDX_SHIFT; // resize vector if necessary int oldSize = offsets.size(); @@ -739,122 +747,98 @@ void BamStandardIndex::SaveLinearOffset(LinearOffsetVector& offsets, for( int i = beginOffset + 1; i <= endOffset; ++i ) { if ( offsets[i] == 0 ) offsets[i] = lastOffset; - } + } } -// initializes index data structure to hold @count references -void BamStandardIndex::SetReferenceCount(const int& count) { - for ( int i = 0; i < count; ++i ) - m_indexData[i].HasAlignments = false; +void BamStandardIndex::SaveLinearOffsetsSummary(const int& refId, const int& numLinearOffsets) { + BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId); + refSummary.NumLinearOffsets = numLinearOffsets; + refSummary.FirstLinearOffsetFilePosition = Tell(); } -bool BamStandardIndex::SkipToFirstReference(void) { - BamStandardIndexData::const_iterator indexBegin = m_indexData.begin(); - return SkipToReference( (*indexBegin).first ); +// seek to position in index file stream +void BamStandardIndex::Seek(const int64_t& position, const int& origin) { + if ( fseek64(Resources.IndexStream, position, origin) != 0 ) + throw BamException("BamStandardIndex::Seek", "could not seek in BAI file"); } -// position file pointer to desired reference begin, return true if skipped OK -bool BamStandardIndex::SkipToReference(const int& refId) { - - // attempt rewind - if ( !Rewind() ) return false; - - // read in number of references - uint32_t numReferences; - size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, m_indexStream); - if ( elementsRead != 1 ) return false; - if ( m_isBigEndian ) SwapEndian_32(numReferences); - - // iterate over reference entries - bool skippedOk = true; - int currentRefId = 0; - while (currentRefId != refId) { - skippedOk &= LoadReference(currentRefId, false); - ++currentRefId; - } - - // return success - return skippedOk; +// change the index caching behavior +void BamStandardIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) { + m_cacheMode = mode; + // do nothing else here ? cache mode will be ignored from now on, most likely } -// write header to new index file -bool BamStandardIndex::WriteHeader(void) { - - size_t elementsWritten = 0; - - // write magic number - elementsWritten += fwrite("BAI\1", sizeof(char), 4, m_indexStream); - - // store offset of beginning of data - m_dataBeginOffset = ftell64(m_indexStream); +void BamStandardIndex::SkipBins(const int& numBins) { + uint32_t binId; + int32_t numAlignmentChunks; + for (int i = 0; i < numBins; ++i) + ReadBinIntoBuffer(binId, numAlignmentChunks); // results & buffer ignored +} - // return success/failure of write - return (elementsWritten == 4); +void BamStandardIndex::SkipLinearOffsets(const int& numLinearOffsets) { + const unsigned int bytesRequested = numLinearOffsets*BamStandardIndex::SIZEOF_LINEAROFFSET; + ReadIntoBuffer(bytesRequested); } -// write index data for all references to new index file -bool BamStandardIndex::WriteAllReferences(void) { +void BamStandardIndex::SortLinearOffsets(BaiLinearOffsetVector& linearOffsets) { + sort( linearOffsets.begin(), linearOffsets.end() ); +} - size_t elementsWritten = 0; +void BamStandardIndex::SummarizeBins(BaiReferenceSummary& refSummary) { - // write number of reference sequences - int32_t numReferenceSeqs = m_indexData.size(); - if ( m_isBigEndian ) SwapEndian_32(numReferenceSeqs); - elementsWritten += fwrite(&numReferenceSeqs, sizeof(numReferenceSeqs), 1, m_indexStream); + // load number of bins + int numBins; + ReadNumBins(numBins); - // iterate over reference sequences - bool refsOk = true; - BamStandardIndexData::const_iterator indexIter = m_indexData.begin(); - BamStandardIndexData::const_iterator indexEnd = m_indexData.end(); - for ( ; indexIter != indexEnd; ++ indexIter ) - refsOk &= WriteReference( (*indexIter).second ); + // store bins summary for this reference + refSummary.NumBins = numBins; + refSummary.FirstBinFilePosition = Tell(); - // return success/failure of write - return ( (elementsWritten == 1) && refsOk ); + // skip this reference's bins + SkipBins(numBins); } -// write index data for bin to new index file -bool BamStandardIndex::WriteBin(const uint32_t& binId, const ChunkVector& chunks) { +void BamStandardIndex::SummarizeIndexFile(void) { - size_t elementsWritten = 0; - - // write BAM bin ID - uint32_t binKey = binId; - if ( m_isBigEndian ) SwapEndian_32(binKey); - elementsWritten += fwrite(&binKey, sizeof(binKey), 1, m_indexStream); + // load number of reference sequences + int numReferences; + ReadNumReferences(numReferences); - // write chunks - bool chunksOk = WriteChunks(chunks); + // initialize file summary data + ReserveForSummary(numReferences); - // return success/failure of write - return ( (elementsWritten == 1) && chunksOk ); + // iterate over reference entries + BaiFileSummary::iterator summaryIter = m_indexFileSummary.begin(); + BaiFileSummary::iterator summaryEnd = m_indexFileSummary.end(); + for ( int i = 0; summaryIter != summaryEnd; ++summaryIter, ++i ) + SummarizeReference(*summaryIter); } -// write index data for bins to new index file -bool BamStandardIndex::WriteBins(const BamBinMap& bins) { +void BamStandardIndex::SummarizeLinearOffsets(BaiReferenceSummary& refSummary) { - size_t elementsWritten = 0; + // load number of linear offsets + int numLinearOffsets; + ReadNumLinearOffsets(numLinearOffsets); - // write number of bins - int32_t binCount = bins.size(); - if ( m_isBigEndian ) SwapEndian_32(binCount); - elementsWritten += fwrite(&binCount, sizeof(binCount), 1, m_indexStream); + // store bin summary data for this reference + refSummary.NumLinearOffsets = numLinearOffsets; + refSummary.FirstLinearOffsetFilePosition = Tell(); - // iterate over bins - bool binsOk = true; - BamBinMap::const_iterator binIter = bins.begin(); - BamBinMap::const_iterator binEnd = bins.end(); - for ( ; binIter != binEnd; ++binIter ) - binsOk &= WriteBin( (*binIter).first, (*binIter).second ); + // skip linear offsets in index file + SkipLinearOffsets(numLinearOffsets); +} - // return success/failure of write - return ( (elementsWritten == 1) && binsOk ); +void BamStandardIndex::SummarizeReference(BaiReferenceSummary& refSummary) { + SummarizeBins(refSummary); + SummarizeLinearOffsets(refSummary); } -// write index data for chunk entry to new index file -bool BamStandardIndex::WriteChunk(const Chunk& chunk) { +// return position of file pointer in index file stream +int64_t BamStandardIndex::Tell(void) const { + return ftell64(Resources.IndexStream); +} - size_t elementsWritten = 0; +void BamStandardIndex::WriteAlignmentChunk(const BaiAlignmentChunk& chunk) { // localize alignment chunk offsets uint64_t start = chunk.Start; @@ -867,63 +851,111 @@ bool BamStandardIndex::WriteChunk(const Chunk& chunk) { } // write to index file - elementsWritten += fwrite(&start, sizeof(start), 1, m_indexStream); - elementsWritten += fwrite(&stop, sizeof(stop), 1, m_indexStream); - - // return success/failure of write - return ( elementsWritten == 2 ); + size_t elementsWritten = 0; + elementsWritten += fwrite(&start, sizeof(start), 1, Resources.IndexStream); + elementsWritten += fwrite(&stop, sizeof(stop), 1, Resources.IndexStream); + if ( elementsWritten != 2 ) + throw BamException("BamStandardIndex::WriteAlignmentChunk", "could not write BAI alignment chunk"); } -// write index data for chunk entry to new index file -bool BamStandardIndex::WriteChunks(const ChunkVector& chunks) { +void BamStandardIndex::WriteAlignmentChunks(BaiAlignmentChunkVector& chunks) { - size_t elementsWritten = 0; + // make sure chunks are merged (simplified) before writing & saving summary + MergeAlignmentChunks(chunks); // write chunks int32_t chunkCount = chunks.size(); if ( m_isBigEndian ) SwapEndian_32(chunkCount); - elementsWritten += fwrite(&chunkCount, sizeof(chunkCount), 1, m_indexStream); + const size_t elementsWritten = fwrite(&chunkCount, sizeof(chunkCount), 1, Resources.IndexStream); + if ( elementsWritten != 1 ) + throw BamException("BamStandardIndex::WriteAlignmentChunks", "could not write BAI chunk count"); // iterate over chunks - bool chunksOk = true; - ChunkVector::const_iterator chunkIter = chunks.begin(); - ChunkVector::const_iterator chunkEnd = chunks.end(); + BaiAlignmentChunkVector::const_iterator chunkIter = chunks.begin(); + BaiAlignmentChunkVector::const_iterator chunkEnd = chunks.end(); for ( ; chunkIter != chunkEnd; ++chunkIter ) - chunksOk &= WriteChunk( (*chunkIter) ); + WriteAlignmentChunk( (*chunkIter) ); +} + +void BamStandardIndex::WriteBin(const uint32_t& binId, BaiAlignmentChunkVector& chunks) { + + // write BAM bin ID + uint32_t binKey = binId; + if ( m_isBigEndian ) SwapEndian_32(binKey); + const size_t elementsWritten = fwrite(&binKey, sizeof(binKey), 1, Resources.IndexStream); + if ( elementsWritten != 1 ) + throw BamException("BamStandardIndex::WriteBin", "could not write bin ID"); - // return success/failure of write - return ( (elementsWritten == 1) && chunksOk ); + // write bin's alignment chunks + WriteAlignmentChunks(chunks); } -// write index data for linear offsets entry to new index file -bool BamStandardIndex::WriteLinearOffsets(const LinearOffsetVector& offsets) { +void BamStandardIndex::WriteBins(const int& refId, BaiBinMap& bins) { + + // write number of bins + int32_t binCount = bins.size(); + if ( m_isBigEndian ) SwapEndian_32(binCount); + const size_t elementsWritten = fwrite(&binCount, sizeof(binCount), 1, Resources.IndexStream); + if ( elementsWritten != 1 ) + throw BamException("BamStandardIndex::WriteBins", "could not write bin count"); + + // save summary for reference's bins + SaveBinsSummary(refId, bins.size()); + + // iterate over bins + BaiBinMap::iterator binIter = bins.begin(); + BaiBinMap::iterator binEnd = bins.end(); + for ( ; binIter != binEnd; ++binIter ) + WriteBin( (*binIter).first, (*binIter).second ); +} + +void BamStandardIndex::WriteHeader(void) { + + size_t elementsWritten = 0; + + // write magic number + elementsWritten += fwrite(BamStandardIndex::BAI_MAGIC, sizeof(char), 4, Resources.IndexStream); + + // write number of reference sequences + int32_t numReferences = m_indexFileSummary.size(); + if ( m_isBigEndian ) SwapEndian_32(numReferences); + elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, Resources.IndexStream); + + if ( elementsWritten != 5 ) + throw BamException("BamStandardIndex::WriteHeader", "could not write BAI header"); +} + +void BamStandardIndex::WriteLinearOffsets(const int& refId, BaiLinearOffsetVector& linearOffsets) { + + // make sure linear offsets are sorted before writing & saving summary + SortLinearOffsets(linearOffsets); size_t elementsWritten = 0; // write number of linear offsets - int32_t offsetCount = offsets.size(); + int32_t offsetCount = linearOffsets.size(); if ( m_isBigEndian ) SwapEndian_32(offsetCount); - elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, m_indexStream); + elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, Resources.IndexStream); + + // save summary for reference's linear offsets + SaveLinearOffsetsSummary(refId, linearOffsets.size()); // iterate over linear offsets - LinearOffsetVector::const_iterator offsetIter = offsets.begin(); - LinearOffsetVector::const_iterator offsetEnd = offsets.end(); + BaiLinearOffsetVector::const_iterator offsetIter = linearOffsets.begin(); + BaiLinearOffsetVector::const_iterator offsetEnd = linearOffsets.end(); for ( ; offsetIter != offsetEnd; ++offsetIter ) { // write linear offset uint64_t linearOffset = (*offsetIter); if ( m_isBigEndian ) SwapEndian_64(linearOffset); - elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, m_indexStream); + elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, Resources.IndexStream); } - // return success/failure of write - return ( elementsWritten == (size_t)(offsetCount + 1) ); + if ( elementsWritten != (linearOffsets.size() + 1) ) + throw BamException("BamStandardIndex::WriteLinearOffsets", "could not write BAI linear offsets"); } -// write index data for a single reference to new index file -bool BamStandardIndex::WriteReference(const ReferenceIndex& refEntry) { - bool refOk = true; - refOk &= WriteBins(refEntry.Bins); - refOk &= WriteLinearOffsets(refEntry.Offsets); - return refOk; +void BamStandardIndex::WriteReferenceEntry(BaiReferenceEntry& refEntry) { + WriteBins(refEntry.ID, refEntry.Bins); + WriteLinearOffsets(refEntry.ID, refEntry.LinearOffsets); }