// ***************************************************************************
// 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 <api/BamAlignment.h>
-#include <api/BamReader.h>
-#include <api/internal/BamReader_p.h>
-#include <api/internal/BamStandardIndex_p.h>
-#include <api/internal/BgzfStream_p.h>
+#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;
#include <cstdlib>
#include <cstring>
#include <algorithm>
-#include <iostream>
-#include <map>
+#include <sstream>
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<uint16_t>& 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<uint16_t>& candidateBins,
+ vector<int64_t>& 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<uint16_t>::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<int64_t>& 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<uint32_t, ChunkVector>::const_iterator binIter = binMap.find(binKey);
- if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) {
-
- // iterate over chunks
- const ChunkVector& chunks = (*binIter).second;
- std::vector<Chunk>::const_iterator chunksIter = chunks.begin();
- std::vector<Chunk>::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<int64_t> 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<int64_t>::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<uint16_t> 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<int64_t> 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<int64_t>::const_iterator OffsetConstIterator;
+ OffsetConstIterator offsetFirst = offsets.begin();
+ OffsetConstIterator offsetIter = offsetFirst;
+ OffsetConstIterator offsetLast = offsets.end();
+ iterator_traits<OffsetConstIterator>::difference_type count = distance(offsetFirst, offsetLast);
+ iterator_traits<OffsetConstIterator>::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<uint32_t, ChunkVector>(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<int32_t, ReferenceIndex>(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<uint32_t, ChunkVector>(saveBin, newChunks));
+ binMap.insert( pair<uint32_t, BaiAlignmentChunkVector>(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();
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;
}
// 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);
}
projects / bamtools.git / blobdiff