1 // ***************************************************************************
2 // BamStandardIndex.cpp (c) 2010 Derek Barnett
3 // Marth Lab, Department of Biology, Boston College
4 // ---------------------------------------------------------------------------
5 // Last modified: 8 October 2011 (DB)
6 // ---------------------------------------------------------------------------
7 // Provides index operations for the standardized BAM index format (".bai")
8 // ***************************************************************************
10 #include <api/BamAlignment.h>
11 #include <api/internal/BamException_p.h>
12 #include <api/internal/BamReader_p.h>
13 #include <api/internal/BamStandardIndex_p.h>
14 using namespace BamTools;
15 using namespace BamTools::Internal;
24 // -----------------------------------
25 // static BamStandardIndex constants
26 // -----------------------------------
28 const int BamStandardIndex::MAX_BIN = 37450; // =(8^6-1)/7+1
29 const int BamStandardIndex::BAM_LIDX_SHIFT = 14;
30 const string BamStandardIndex::BAI_EXTENSION = ".bai";
31 const char* const BamStandardIndex::BAI_MAGIC = "BAI\1";
32 const int BamStandardIndex::SIZEOF_ALIGNMENTCHUNK = sizeof(uint64_t)*2;
33 const int BamStandardIndex::SIZEOF_BINCORE = sizeof(uint32_t) + sizeof(int32_t);
34 const int BamStandardIndex::SIZEOF_LINEAROFFSET = sizeof(uint64_t);
36 // ----------------------------
37 // RaiiWrapper implementation
38 // ----------------------------
40 BamStandardIndex::RaiiWrapper::RaiiWrapper(void)
45 BamStandardIndex::RaiiWrapper::~RaiiWrapper(void) {
58 // ---------------------------------
59 // BamStandardIndex implementation
60 // ---------------------------------
63 BamStandardIndex::BamStandardIndex(Internal::BamReaderPrivate* reader)
65 , m_cacheMode(BamIndex::LimitedIndexCaching)
68 m_isBigEndian = BamTools::SystemIsBigEndian();
72 BamStandardIndex::~BamStandardIndex(void) {
76 void BamStandardIndex::AdjustRegion(const BamRegion& region, uint32_t& begin, uint32_t& end) {
78 // retrieve references from reader
79 const RefVector& references = m_reader->GetReferenceData();
81 // LeftPosition cannot be greater than or equal to reference length
82 if ( region.LeftPosition >= references.at(region.LeftRefID).RefLength )
83 throw BamException("BamStandardIndex::AdjustRegion", "invalid region requested");
86 begin = (unsigned int)region.LeftPosition;
88 // if right bound specified AND left&right bounds are on same reference
89 // OK to use right bound position as region 'end'
90 if ( region.isRightBoundSpecified() && ( region.LeftRefID == region.RightRefID ) )
91 end = (unsigned int)region.RightPosition;
93 // otherwise, set region 'end' to last reference base
94 else end = (unsigned int)references.at(region.LeftRefID).RefLength;
98 void BamStandardIndex::CalculateCandidateBins(const uint32_t& begin,
100 set<uint16_t>& candidateBins)
102 // initialize list, bin '0' is always a valid bin
103 candidateBins.insert(0);
105 // get rest of bins that contain this region
107 for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { candidateBins.insert(k); }
108 for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { candidateBins.insert(k); }
109 for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { candidateBins.insert(k); }
110 for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { candidateBins.insert(k); }
111 for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { candidateBins.insert(k); }
114 void BamStandardIndex::CalculateCandidateOffsets(const BaiReferenceSummary& refSummary,
115 const uint64_t& minOffset,
116 set<uint16_t>& candidateBins,
117 vector<int64_t>& offsets)
120 Seek(refSummary.FirstBinFilePosition, SEEK_SET);
122 // iterate over reference bins
124 int32_t numAlignmentChunks;
125 set<uint16_t>::iterator candidateBinIter;
126 for ( int i = 0; i < refSummary.NumBins; ++i ) {
128 // read bin contents (if successful, alignment chunks are now in m_buffer)
129 ReadBinIntoBuffer(binId, numAlignmentChunks);
131 // see if bin is a 'candidate bin'
132 candidateBinIter = candidateBins.find(binId);
134 // if not, move on to next bin
135 if ( candidateBinIter == candidateBins.end() )
138 // otherwise, check bin's contents against for overlap
145 // iterate over alignment chunks
146 for ( int j = 0; j < numAlignmentChunks; ++j ) {
148 // read chunk start & stop from buffer
149 memcpy((char*)&chunkStart, Resources.Buffer+offset, sizeof(uint64_t));
150 offset += sizeof(uint64_t);
151 memcpy((char*)&chunkStop, Resources.Buffer+offset, sizeof(uint64_t));
152 offset += sizeof(uint64_t);
154 // swap endian-ness if necessary
155 if ( m_isBigEndian ) {
156 SwapEndian_64(chunkStart);
157 SwapEndian_64(chunkStop);
160 // store alignment chunk's start offset
161 // if its stop offset is larger than our 'minOffset'
162 if ( chunkStop >= minOffset )
163 offsets.push_back(chunkStart);
166 // 'pop' bin ID from candidate bins set
167 candidateBins.erase(candidateBinIter);
169 // quit if no more candidates
170 if ( candidateBins.empty() )
176 uint64_t BamStandardIndex::CalculateMinOffset(const BaiReferenceSummary& refSummary,
177 const uint32_t& begin)
179 // if no linear offsets exist, return 0
180 if ( refSummary.NumLinearOffsets == 0 )
183 // if 'begin' starts beyond last linear offset, use the last linear offset as minimum
184 // else use the offset corresponding to the requested start position
185 const int shiftedBegin = begin>>BamStandardIndex::BAM_LIDX_SHIFT;
186 if ( shiftedBegin >= refSummary.NumLinearOffsets )
187 return LookupLinearOffset( refSummary, refSummary.NumLinearOffsets-1 );
189 return LookupLinearOffset( refSummary, shiftedBegin );
192 void BamStandardIndex::CheckBufferSize(char*& buffer,
193 unsigned int& bufferLength,
194 const unsigned int& requestedBytes)
197 if ( requestedBytes > bufferLength ) {
198 bufferLength = requestedBytes + 10;
200 buffer = new char[bufferLength];
202 } catch ( std::bad_alloc& ) {
204 s << "out of memory when allocating " << requestedBytes << " bytes";
205 throw BamException("BamStandardIndex::CheckBufferSize", s.str());
209 void BamStandardIndex::CheckBufferSize(unsigned char*& buffer,
210 unsigned int& bufferLength,
211 const unsigned int& requestedBytes)
214 if ( requestedBytes > bufferLength ) {
215 bufferLength = requestedBytes + 10;
217 buffer = new unsigned char[bufferLength];
219 } catch ( std::bad_alloc& ) {
221 s << "out of memory when allocating " << requestedBytes << " bytes";
222 throw BamException("BamStandardIndex::CheckBufferSize", s.str());
226 void BamStandardIndex::CheckMagicNumber(void) {
228 // check 'magic number' to see if file is BAI index
230 const size_t elementsRead = fread(magic, sizeof(char), 4, Resources.IndexStream);
231 if ( elementsRead != 4 )
232 throw BamException("BamStandardIndex::CheckMagicNumber", "could not read BAI magic number");
234 // compare to expected value
235 if ( strncmp(magic, BamStandardIndex::BAI_MAGIC, 4) != 0 )
236 throw BamException("BamStandardIndex::CheckMagicNumber", "invalid BAI magic number");
239 void BamStandardIndex::ClearReferenceEntry(BaiReferenceEntry& refEntry) {
241 refEntry.Bins.clear();
242 refEntry.LinearOffsets.clear();
245 void BamStandardIndex::CloseFile(void) {
248 if ( IsFileOpen() ) {
249 fclose(Resources.IndexStream);
250 Resources.IndexStream = 0;
253 // clear index file summary data
254 m_indexFileSummary.clear();
256 // clean up I/O buffer
257 delete[] Resources.Buffer;
258 Resources.Buffer = 0;
262 // builds index from associated BAM file & writes out to index file
263 bool BamStandardIndex::Create(void) {
265 // skip if BamReader is invalid or not open
266 if ( m_reader == 0 || !m_reader->IsOpen() ) {
267 SetErrorString("BamStandardIndex::Create", "could not create index: reader is not open");
272 if ( !m_reader->Rewind() ) {
273 const string readerError = m_reader->GetErrorString();
274 const string message = "could not create index: \n\t" + readerError;
275 SetErrorString("BamStandardIndex::Create", message);
281 // open new index file (read & write)
282 string indexFilename = m_reader->Filename() + Extension();
283 OpenFile(indexFilename, "w+b");
285 // initialize BaiFileSummary with number of references
286 const int& numReferences = m_reader->GetReferenceCount();
287 ReserveForSummary(numReferences);
289 // initialize output file
292 // set up bin, ID, offset, & coordinate markers
293 const uint32_t defaultValue = 0xffffffffu;
294 uint32_t currentBin = defaultValue;
295 uint32_t lastBin = defaultValue;
296 int32_t currentRefID = defaultValue;
297 int32_t lastRefID = defaultValue;
298 uint64_t currentOffset = (uint64_t)m_reader->Tell();
299 uint64_t lastOffset = currentOffset;
300 int32_t lastPosition = defaultValue;
302 // iterate through alignments in BAM file
304 BaiReferenceEntry refEntry;
305 while ( m_reader->LoadNextAlignment(al) ) {
307 // changed to new reference
308 if ( lastRefID != al.RefID ) {
310 // if not first reference, save previous reference data
311 if ( lastRefID != (int32_t)defaultValue ) {
313 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
314 WriteReferenceEntry(refEntry);
315 ClearReferenceEntry(refEntry);
317 // write any empty references between (but *NOT* including) lastRefID & al.RefID
318 for ( int i = lastRefID+1; i < al.RefID; ++i ) {
319 BaiReferenceEntry emptyEntry(i);
320 WriteReferenceEntry(emptyEntry);
323 // update bin markers
324 currentOffset = lastOffset;
327 currentRefID = al.RefID;
330 // otherwise, this is first pass
331 // be sure to write any empty references up to (but *NOT* including) current RefID
333 for ( int i = 0; i < al.RefID; ++i ) {
334 BaiReferenceEntry emptyEntry(i);
335 WriteReferenceEntry(emptyEntry);
339 // update reference markers
340 refEntry.ID = al.RefID;
341 lastRefID = al.RefID;
342 lastBin = defaultValue;
345 // if lastPosition greater than current alignment position - file not sorted properly
346 else if ( lastPosition > al.Position ) {
348 s << "BAM file is not properly sorted by coordinate" << endl
349 << "Current alignment position: " << al.Position
350 << " < previous alignment position: " << lastPosition
351 << " on reference ID: " << al.RefID << endl;
352 SetErrorString("BamStandardIndex::Create", s.str());
356 // if alignment's ref ID is valid & its bin is not a 'leaf'
357 if ( (al.RefID >= 0) && (al.Bin < 4681) )
358 SaveLinearOffsetEntry(refEntry.LinearOffsets, al.Position, al.GetEndPosition(), lastOffset);
360 // changed to new BAI bin
361 if ( al.Bin != lastBin ) {
363 // if not first bin on reference, save previous bin data
364 if ( currentBin != defaultValue )
365 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
368 currentOffset = lastOffset;
371 currentRefID = al.RefID;
373 // if invalid RefID, break out
374 if ( currentRefID < 0 )
378 // make sure that current file pointer is beyond lastOffset
379 if ( m_reader->Tell() <= (int64_t)lastOffset ) {
380 SetErrorString("BamStandardIndex::Create", "calculating offsets failed");
384 // update lastOffset & lastPosition
385 lastOffset = m_reader->Tell();
386 lastPosition = al.Position;
389 // after finishing alignments, if any data was read, check:
390 if ( currentRefID >= 0 ) {
392 // store last alignment chunk to its bin, then write last reference entry with data
393 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
394 WriteReferenceEntry(refEntry);
396 // then write any empty references remaining at end of file
397 for ( int i = currentRefID+1; i < numReferences; ++i ) {
398 BaiReferenceEntry emptyEntry(i);
399 WriteReferenceEntry(emptyEntry);
403 } catch ( BamException& e) {
404 m_errorString = e.what();
409 if ( !m_reader->Rewind() ) {
410 const string readerError = m_reader->GetErrorString();
411 const string message = "could not create index: \n\t" + readerError;
412 SetErrorString("BamStandardIndex::Create", message);
420 // returns format's file extension
421 const string BamStandardIndex::Extension(void) {
422 return BamStandardIndex::BAI_EXTENSION;
425 void BamStandardIndex::GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion) {
427 // cannot calculate offsets if unknown/invalid reference ID requested
428 if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() )
429 throw BamException("BamStandardIndex::GetOffset", "invalid reference ID requested");
431 // retrieve index summary for left bound reference
432 const BaiReferenceSummary& refSummary = m_indexFileSummary.at(region.LeftRefID);
434 // set up region boundaries based on actual BamReader data
437 AdjustRegion(region, begin, end);
439 // retrieve all candidate bin IDs for region
440 set<uint16_t> candidateBins;
441 CalculateCandidateBins(begin, end, candidateBins);
443 // use reference's linear offsets to calculate the minimum offset
444 // that must be considered to find overlap
445 const uint64_t& minOffset = CalculateMinOffset(refSummary, begin);
447 // attempt to use reference summary, minOffset, & candidateBins to calculate offsets
448 // no data should not be error, just bail
449 vector<int64_t> offsets;
450 CalculateCandidateOffsets(refSummary, minOffset, candidateBins, offsets);
451 if ( offsets.empty() )
454 // ensure that offsets are sorted before processing
455 sort( offsets.begin(), offsets.end() );
457 // binary search for an overlapping block (may not be first one though)
459 typedef vector<int64_t>::const_iterator OffsetConstIterator;
460 OffsetConstIterator offsetFirst = offsets.begin();
461 OffsetConstIterator offsetIter = offsetFirst;
462 OffsetConstIterator offsetLast = offsets.end();
463 iterator_traits<OffsetConstIterator>::difference_type count = distance(offsetFirst, offsetLast);
464 iterator_traits<OffsetConstIterator>::difference_type step;
465 while ( count > 0 ) {
466 offsetIter = offsetFirst;
468 advance(offsetIter, step);
470 // attempt seek to candidate offset
471 const int64_t& candidateOffset = (*offsetIter);
472 if ( !m_reader->Seek(candidateOffset) ) {
473 const string readerError = m_reader->GetErrorString();
474 const string message = "could not seek in BAM file: \n\t" + readerError;
475 throw BamException("BamToolsIndex::GetOffset", message);
478 // load first available alignment, setting flag to true if data exists
479 *hasAlignmentsInRegion = m_reader->LoadNextAlignment(al);
481 // check alignment against region
482 if ( al.GetEndPosition() <= region.LeftPosition ) {
483 offsetFirst = ++offsetIter;
488 // step back to the offset before the 'current offset' (to make sure we cover overlaps)
489 if ( offsetIter != offsets.begin() )
491 offset = (*offsetIter);
494 // returns whether reference has alignments or no
495 bool BamStandardIndex::HasAlignments(const int& referenceID) const {
496 if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() )
498 const BaiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID);
499 return ( refSummary.NumBins > 0 );
502 bool BamStandardIndex::IsFileOpen(void) const {
503 return ( Resources.IndexStream != 0 );
506 // attempts to use index data to jump to @region, returns success/fail
507 // a "successful" jump indicates no error, but not whether this region has data
508 // * thus, the method sets a flag to indicate whether there are alignments
509 // available after the jump position
510 bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
513 *hasAlignmentsInRegion = false;
515 // skip if invalid reader or not open
516 if ( m_reader == 0 || !m_reader->IsOpen() ) {
517 SetErrorString("BamStandardIndex::Jump", "could not jump: reader is not open");
521 // calculate nearest offset to jump to
524 GetOffset(region, offset, hasAlignmentsInRegion);
525 } catch ( BamException& e ) {
526 m_errorString = e.what();
530 // if region has alignments, return success/fail of seeking there
531 if ( *hasAlignmentsInRegion )
532 return m_reader->Seek(offset);
534 // otherwise, simply return true (but hasAlignmentsInRegion flag has been set to false)
535 // (this is OK, BamReader will check this flag before trying to load data)
539 // loads existing data from file into memory
540 bool BamStandardIndex::Load(const std::string& filename) {
544 // attempt to open file (read-only)
545 OpenFile(filename, "rb");
550 // load in-memory summary of index data
551 SummarizeIndexFile();
556 } catch ( BamException& e ) {
557 m_errorString = e.what();
562 uint64_t BamStandardIndex::LookupLinearOffset(const BaiReferenceSummary& refSummary, const int& index) {
564 // attempt seek to proper index file position
565 const int64_t linearOffsetFilePosition = (int64_t)refSummary.FirstLinearOffsetFilePosition +
566 index*BamStandardIndex::SIZEOF_LINEAROFFSET;
567 Seek(linearOffsetFilePosition, SEEK_SET);
569 // read linear offset from BAI file
570 uint64_t linearOffset;
571 ReadLinearOffset(linearOffset);
575 void BamStandardIndex::MergeAlignmentChunks(BaiAlignmentChunkVector& chunks) {
577 // skip if chunks are empty, nothing to merge
578 if ( chunks.empty() )
581 // set up merged alignment chunk container
582 BaiAlignmentChunkVector mergedChunks;
583 mergedChunks.push_back( chunks[0] );
585 // iterate over chunks
587 BaiAlignmentChunkVector::iterator chunkIter = chunks.begin();
588 BaiAlignmentChunkVector::iterator chunkEnd = chunks.end();
589 for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {
591 // get 'currentMergeChunk' based on numeric index
592 BaiAlignmentChunk& currentMergeChunk = mergedChunks[i];
594 // get sourceChunk based on source vector iterator
595 BaiAlignmentChunk& sourceChunk = (*chunkIter);
597 // if currentMergeChunk ends where sourceChunk starts, then merge the two
598 if ( currentMergeChunk.Stop>>16 == sourceChunk.Start>>16 )
599 currentMergeChunk.Stop = sourceChunk.Stop;
603 // append sourceChunk after currentMergeChunk
604 mergedChunks.push_back(sourceChunk);
606 // update i, so the next iteration will consider the
607 // recently-appended sourceChunk as new mergeChunk candidate
612 // saved newly-merged chunks into (parameter) chunks
613 chunks = mergedChunks;
616 void BamStandardIndex::OpenFile(const std::string& filename, const char* mode) {
618 // make sure any previous index file is closed
621 // attempt to open file
622 Resources.IndexStream = fopen(filename.c_str(), mode);
623 if ( !IsFileOpen() ) {
624 const string message = string("could not open file: ") + filename;
625 throw BamException("BamStandardIndex::OpenFile", message);
629 void BamStandardIndex::ReadBinID(uint32_t& binId) {
630 const size_t elementsRead = fread(&binId, sizeof(binId), 1, Resources.IndexStream);
631 if ( m_isBigEndian ) SwapEndian_32(binId);
632 if ( elementsRead != 1 )
633 throw BamException("BamStandardIndex::ReadBinID", "could not read BAI bin ID");
636 void BamStandardIndex::ReadBinIntoBuffer(uint32_t& binId, int32_t& numAlignmentChunks) {
640 ReadNumAlignmentChunks(numAlignmentChunks);
643 const unsigned int bytesRequested = numAlignmentChunks*BamStandardIndex::SIZEOF_ALIGNMENTCHUNK;
644 ReadIntoBuffer(bytesRequested);
647 void BamStandardIndex::ReadIntoBuffer(const unsigned int& bytesRequested) {
649 // ensure that our buffer is big enough for request
650 BamStandardIndex::CheckBufferSize(Resources.Buffer, m_bufferLength, bytesRequested);
652 // read from BAI file stream
653 const size_t bytesRead = fread( Resources.Buffer, sizeof(char), bytesRequested, Resources.IndexStream );
654 if ( bytesRead != (size_t)bytesRequested ) {
656 s << "expected to read: " << bytesRequested << " bytes, "
657 << "but instead read: " << bytesRead;
658 throw BamException("BamStandardIndex::ReadIntoBuffer", s.str());
662 void BamStandardIndex::ReadLinearOffset(uint64_t& linearOffset) {
663 const size_t elementsRead = fread(&linearOffset, sizeof(linearOffset), 1, Resources.IndexStream);
664 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
665 if ( elementsRead != 1 )
666 throw BamException("BamStandardIndex::ReadLinearOffset", "could not read BAI linear offset");
669 void BamStandardIndex::ReadNumAlignmentChunks(int& numAlignmentChunks) {
670 const size_t elementsRead = fread(&numAlignmentChunks, sizeof(numAlignmentChunks), 1, Resources.IndexStream);
671 if ( m_isBigEndian ) SwapEndian_32(numAlignmentChunks);
672 if ( elementsRead != 1 )
673 throw BamException("BamStandardIndex::ReadNumAlignmentChunks", "could not read BAI chunk count");
676 void BamStandardIndex::ReadNumBins(int& numBins) {
677 const size_t elementsRead = fread(&numBins, sizeof(numBins), 1, Resources.IndexStream);
678 if ( m_isBigEndian ) SwapEndian_32(numBins);
679 if ( elementsRead != 1 )
680 throw BamException("BamStandardIndex::ReadNumBins", "could not read BAI bin count");
683 void BamStandardIndex::ReadNumLinearOffsets(int& numLinearOffsets) {
684 const size_t elementsRead = fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, Resources.IndexStream);
685 if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets);
686 if ( elementsRead != 1 )
687 throw BamException("BamStandardIndex::ReadNumAlignmentChunks", "could not read BAI linear offset count");
690 void BamStandardIndex::ReadNumReferences(int& numReferences) {
691 const size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, Resources.IndexStream);
692 if ( m_isBigEndian ) SwapEndian_32(numReferences);
693 if ( elementsRead != 1 )
694 throw BamException("BamStandardIndex::ReadNumReferences", "could not read reference count");
697 void BamStandardIndex::ReserveForSummary(const int& numReferences) {
698 m_indexFileSummary.clear();
699 m_indexFileSummary.assign( numReferences, BaiReferenceSummary() );
702 void BamStandardIndex::SaveAlignmentChunkToBin(BaiBinMap& binMap,
703 const uint32_t& currentBin,
704 const uint64_t& currentOffset,
705 const uint64_t& lastOffset)
707 // create new alignment chunk
708 BaiAlignmentChunk newChunk(currentOffset, lastOffset);
710 // if no entry exists yet for this bin, create one and store alignment chunk
711 BaiBinMap::iterator binIter = binMap.find(currentBin);
712 if ( binIter == binMap.end() ) {
713 BaiAlignmentChunkVector newChunks;
714 newChunks.push_back(newChunk);
715 binMap.insert( pair<uint32_t, BaiAlignmentChunkVector>(currentBin, newChunks));
718 // otherwise, just append alignment chunk
720 BaiAlignmentChunkVector& binChunks = (*binIter).second;
721 binChunks.push_back( newChunk );
725 void BamStandardIndex::SaveBinsSummary(const int& refId, const int& numBins) {
726 BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
727 refSummary.NumBins = numBins;
728 refSummary.FirstBinFilePosition = Tell();
731 void BamStandardIndex::SaveLinearOffsetEntry(BaiLinearOffsetVector& offsets,
732 const int& alignmentStartPosition,
733 const int& alignmentStopPosition,
734 const uint64_t& lastOffset)
736 // get converted offsets
737 const int beginOffset = alignmentStartPosition >> BamStandardIndex::BAM_LIDX_SHIFT;
738 const int endOffset = (alignmentStopPosition - 1) >> BamStandardIndex::BAM_LIDX_SHIFT;
740 // resize vector if necessary
741 int oldSize = offsets.size();
742 int newSize = endOffset + 1;
743 if ( oldSize < newSize )
744 offsets.resize(newSize, 0);
747 for( int i = beginOffset + 1; i <= endOffset; ++i ) {
748 if ( offsets[i] == 0 )
749 offsets[i] = lastOffset;
753 void BamStandardIndex::SaveLinearOffsetsSummary(const int& refId, const int& numLinearOffsets) {
754 BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
755 refSummary.NumLinearOffsets = numLinearOffsets;
756 refSummary.FirstLinearOffsetFilePosition = Tell();
759 // seek to position in index file stream
760 void BamStandardIndex::Seek(const int64_t& position, const int& origin) {
761 if ( fseek64(Resources.IndexStream, position, origin) != 0 )
762 throw BamException("BamStandardIndex::Seek", "could not seek in BAI file");
765 // change the index caching behavior
766 void BamStandardIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) {
768 // do nothing else here ? cache mode will be ignored from now on, most likely
771 void BamStandardIndex::SkipBins(const int& numBins) {
773 int32_t numAlignmentChunks;
774 for (int i = 0; i < numBins; ++i)
775 ReadBinIntoBuffer(binId, numAlignmentChunks); // results & buffer ignored
778 void BamStandardIndex::SkipLinearOffsets(const int& numLinearOffsets) {
779 const unsigned int bytesRequested = numLinearOffsets*BamStandardIndex::SIZEOF_LINEAROFFSET;
780 ReadIntoBuffer(bytesRequested);
783 void BamStandardIndex::SortLinearOffsets(BaiLinearOffsetVector& linearOffsets) {
784 sort( linearOffsets.begin(), linearOffsets.end() );
787 void BamStandardIndex::SummarizeBins(BaiReferenceSummary& refSummary) {
789 // load number of bins
791 ReadNumBins(numBins);
793 // store bins summary for this reference
794 refSummary.NumBins = numBins;
795 refSummary.FirstBinFilePosition = Tell();
797 // skip this reference's bins
801 void BamStandardIndex::SummarizeIndexFile(void) {
803 // load number of reference sequences
805 ReadNumReferences(numReferences);
807 // initialize file summary data
808 ReserveForSummary(numReferences);
810 // iterate over reference entries
811 BaiFileSummary::iterator summaryIter = m_indexFileSummary.begin();
812 BaiFileSummary::iterator summaryEnd = m_indexFileSummary.end();
813 for ( int i = 0; summaryIter != summaryEnd; ++summaryIter, ++i )
814 SummarizeReference(*summaryIter);
817 void BamStandardIndex::SummarizeLinearOffsets(BaiReferenceSummary& refSummary) {
819 // load number of linear offsets
820 int numLinearOffsets;
821 ReadNumLinearOffsets(numLinearOffsets);
823 // store bin summary data for this reference
824 refSummary.NumLinearOffsets = numLinearOffsets;
825 refSummary.FirstLinearOffsetFilePosition = Tell();
827 // skip linear offsets in index file
828 SkipLinearOffsets(numLinearOffsets);
831 void BamStandardIndex::SummarizeReference(BaiReferenceSummary& refSummary) {
832 SummarizeBins(refSummary);
833 SummarizeLinearOffsets(refSummary);
836 // return position of file pointer in index file stream
837 int64_t BamStandardIndex::Tell(void) const {
838 return ftell64(Resources.IndexStream);
841 void BamStandardIndex::WriteAlignmentChunk(const BaiAlignmentChunk& chunk) {
843 // localize alignment chunk offsets
844 uint64_t start = chunk.Start;
845 uint64_t stop = chunk.Stop;
847 // swap endian-ness if necessary
848 if ( m_isBigEndian ) {
849 SwapEndian_64(start);
853 // write to index file
854 size_t elementsWritten = 0;
855 elementsWritten += fwrite(&start, sizeof(start), 1, Resources.IndexStream);
856 elementsWritten += fwrite(&stop, sizeof(stop), 1, Resources.IndexStream);
857 if ( elementsWritten != 2 )
858 throw BamException("BamStandardIndex::WriteAlignmentChunk", "could not write BAI alignment chunk");
861 void BamStandardIndex::WriteAlignmentChunks(BaiAlignmentChunkVector& chunks) {
863 // make sure chunks are merged (simplified) before writing & saving summary
864 MergeAlignmentChunks(chunks);
867 int32_t chunkCount = chunks.size();
868 if ( m_isBigEndian ) SwapEndian_32(chunkCount);
869 const size_t elementsWritten = fwrite(&chunkCount, sizeof(chunkCount), 1, Resources.IndexStream);
870 if ( elementsWritten != 1 )
871 throw BamException("BamStandardIndex::WriteAlignmentChunks", "could not write BAI chunk count");
873 // iterate over chunks
874 BaiAlignmentChunkVector::const_iterator chunkIter = chunks.begin();
875 BaiAlignmentChunkVector::const_iterator chunkEnd = chunks.end();
876 for ( ; chunkIter != chunkEnd; ++chunkIter )
877 WriteAlignmentChunk( (*chunkIter) );
880 void BamStandardIndex::WriteBin(const uint32_t& binId, BaiAlignmentChunkVector& chunks) {
883 uint32_t binKey = binId;
884 if ( m_isBigEndian ) SwapEndian_32(binKey);
885 const size_t elementsWritten = fwrite(&binKey, sizeof(binKey), 1, Resources.IndexStream);
886 if ( elementsWritten != 1 )
887 throw BamException("BamStandardIndex::WriteBin", "could not write bin ID");
889 // write bin's alignment chunks
890 WriteAlignmentChunks(chunks);
893 void BamStandardIndex::WriteBins(const int& refId, BaiBinMap& bins) {
895 // write number of bins
896 int32_t binCount = bins.size();
897 if ( m_isBigEndian ) SwapEndian_32(binCount);
898 const size_t elementsWritten = fwrite(&binCount, sizeof(binCount), 1, Resources.IndexStream);
899 if ( elementsWritten != 1 )
900 throw BamException("BamStandardIndex::WriteBins", "could not write bin count");
902 // save summary for reference's bins
903 SaveBinsSummary(refId, bins.size());
906 BaiBinMap::iterator binIter = bins.begin();
907 BaiBinMap::iterator binEnd = bins.end();
908 for ( ; binIter != binEnd; ++binIter )
909 WriteBin( (*binIter).first, (*binIter).second );
912 void BamStandardIndex::WriteHeader(void) {
914 size_t elementsWritten = 0;
916 // write magic number
917 elementsWritten += fwrite(BamStandardIndex::BAI_MAGIC, sizeof(char), 4, Resources.IndexStream);
919 // write number of reference sequences
920 int32_t numReferences = m_indexFileSummary.size();
921 if ( m_isBigEndian ) SwapEndian_32(numReferences);
922 elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, Resources.IndexStream);
924 if ( elementsWritten != 5 )
925 throw BamException("BamStandardIndex::WriteHeader", "could not write BAI header");
928 void BamStandardIndex::WriteLinearOffsets(const int& refId, BaiLinearOffsetVector& linearOffsets) {
930 // make sure linear offsets are sorted before writing & saving summary
931 SortLinearOffsets(linearOffsets);
933 size_t elementsWritten = 0;
935 // write number of linear offsets
936 int32_t offsetCount = linearOffsets.size();
937 if ( m_isBigEndian ) SwapEndian_32(offsetCount);
938 elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, Resources.IndexStream);
940 // save summary for reference's linear offsets
941 SaveLinearOffsetsSummary(refId, linearOffsets.size());
943 // iterate over linear offsets
944 BaiLinearOffsetVector::const_iterator offsetIter = linearOffsets.begin();
945 BaiLinearOffsetVector::const_iterator offsetEnd = linearOffsets.end();
946 for ( ; offsetIter != offsetEnd; ++offsetIter ) {
948 // write linear offset
949 uint64_t linearOffset = (*offsetIter);
950 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
951 elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, Resources.IndexStream);
954 if ( elementsWritten != (linearOffsets.size() + 1) )
955 throw BamException("BamStandardIndex::WriteLinearOffsets", "could not write BAI linear offsets");
958 void BamStandardIndex::WriteReferenceEntry(BaiReferenceEntry& refEntry) {
959 WriteBins(refEntry.ID, refEntry.Bins);
960 WriteLinearOffsets(refEntry.ID, refEntry.LinearOffsets);