1 // ***************************************************************************
2 // BamStandardIndex.cpp (c) 2010 Derek Barnett
3 // Marth Lab, Department of Biology, Boston College
4 // All rights reserved.
5 // ---------------------------------------------------------------------------
6 // Last modified: 27 April 2011 (DB)
7 // ---------------------------------------------------------------------------
8 // Provides index operations for the standardized BAM index format (".bai")
9 // ***************************************************************************
11 #include <api/BamAlignment.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 // static BamStandardIndex constants
25 const int BamStandardIndex::MAX_BIN = 37450; // =(8^6-1)/7+1
26 const int BamStandardIndex::BAM_LIDX_SHIFT = 14;
27 const string BamStandardIndex::BAI_EXTENSION = ".bai";
28 const char* const BamStandardIndex::BAI_MAGIC = "BAI\1";
29 const int BamStandardIndex::SIZEOF_ALIGNMENTCHUNK = sizeof(uint64_t)*2;
30 const int BamStandardIndex::SIZEOF_BINCORE = sizeof(uint32_t) + sizeof(int32_t);
31 const int BamStandardIndex::SIZEOF_LINEAROFFSET = sizeof(uint64_t);
34 BamStandardIndex::BamStandardIndex(Internal::BamReaderPrivate* reader)
37 , m_cacheMode(BamIndex::LimitedIndexCaching)
41 m_isBigEndian = BamTools::SystemIsBigEndian();
45 BamStandardIndex::~BamStandardIndex(void) {
49 bool BamStandardIndex::AdjustRegion(const BamRegion& region, uint32_t& begin, uint32_t& end) {
51 // retrieve references from reader
52 const RefVector& references = m_reader->GetReferenceData();
54 // make sure left-bound position is valid
55 if ( region.LeftPosition > references.at(region.LeftRefID).RefLength )
59 begin = (unsigned int)region.LeftPosition;
61 // if right bound specified AND left&right bounds are on same reference
62 // OK to use right bound position as region 'end'
63 if ( region.isRightBoundSpecified() && ( region.LeftRefID == region.RightRefID ) )
64 end = (unsigned int)region.RightPosition;
66 // otherwise, set region 'end' to last reference base
67 else end = (unsigned int)references.at(region.LeftRefID).RefLength - 1;
73 void BamStandardIndex::CalculateCandidateBins(const uint32_t& begin,
75 set<uint16_t>& candidateBins)
77 // initialize list, bin '0' is always a valid bin
78 candidateBins.insert(0);
80 // get rest of bins that contain this region
82 for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { candidateBins.insert(k); }
83 for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { candidateBins.insert(k); }
84 for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { candidateBins.insert(k); }
85 for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { candidateBins.insert(k); }
86 for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { candidateBins.insert(k); }
89 bool BamStandardIndex::CalculateCandidateOffsets(const BaiReferenceSummary& refSummary,
90 const uint64_t& minOffset,
91 set<uint16_t>& candidateBins,
92 vector<int64_t>& offsets)
94 // attempt seek to first bin
95 if ( !Seek(refSummary.FirstBinFilePosition, SEEK_SET) )
98 // iterate over reference bins
100 int32_t numAlignmentChunks;
101 set<uint16_t>::iterator candidateBinIter;
102 for ( int i = 0; i < refSummary.NumBins; ++i ) {
104 // read bin contents (if successful, alignment chunks are now in m_buffer)
105 if ( !ReadBinIntoBuffer(binId, numAlignmentChunks) )
108 // see if bin is a 'candidate bin'
109 candidateBinIter = candidateBins.find(binId);
111 // if not, move on to next bin
112 if ( candidateBinIter == candidateBins.end() )
115 // otherwise, check bin's contents against for overlap
118 unsigned int offset = 0;
122 // iterate over alignment chunks
123 for (int j = 0; j < numAlignmentChunks; ++j ) {
125 // read chunk start & stop from buffer
126 memcpy((char*)&chunkStart, m_buffer+offset, sizeof(uint64_t));
127 offset += sizeof(uint64_t);
128 memcpy((char*)&chunkStop, m_buffer, sizeof(uint64_t));
129 offset += sizeof(uint64_t);
131 // swap endian-ness if necessary
132 if ( m_isBigEndian ) {
133 SwapEndian_64(chunkStart);
134 SwapEndian_64(chunkStop);
137 // store alignment chunk's start offset
138 // if its stop offset is larger than our 'minOffset'
139 if ( chunkStop > minOffset )
140 offsets.push_back(chunkStart);
143 // 'pop' bin ID from candidate bins set
144 candidateBins.erase(candidateBinIter);
146 // quit if no more candidates
147 if ( candidateBins.empty() )
152 // return success/failure on calculating at least 1 offset
153 return ( !offsets.empty() );
156 uint64_t BamStandardIndex::CalculateMinOffset(const BaiReferenceSummary& refSummary,
157 const uint32_t& begin)
159 // if no linear offsets exist, return 0
160 if ( refSummary.NumLinearOffsets == 0 )
163 // if 'begin' starts beyond last linear offset, use the last linear offset as minimum
164 // else use the offset corresponding to the requested start position
165 const int shiftedBegin = begin>>BamStandardIndex::BAM_LIDX_SHIFT;
166 if ( shiftedBegin >= refSummary.NumLinearOffsets )
167 return LookupLinearOffset( refSummary, refSummary.NumLinearOffsets-1 );
169 return LookupLinearOffset( refSummary, shiftedBegin );
172 void BamStandardIndex::CheckBufferSize(char*& buffer,
173 unsigned int& bufferLength,
174 const unsigned int& requestedBytes)
177 if ( requestedBytes > bufferLength ) {
178 bufferLength = requestedBytes + 10;
180 buffer = new char[bufferLength];
182 } catch ( std::bad_alloc ) {
183 cerr << "BamStandardIndex ERROR: out of memory when allocating "
184 << requestedBytes << " byes" << endl;
189 void BamStandardIndex::CheckBufferSize(unsigned char*& buffer,
190 unsigned int& bufferLength,
191 const unsigned int& requestedBytes)
194 if ( requestedBytes > bufferLength ) {
195 bufferLength = requestedBytes + 10;
197 buffer = new unsigned char[bufferLength];
199 } catch ( std::bad_alloc ) {
200 cerr << "BamStandardIndex ERROR: out of memory when allocating "
201 << requestedBytes << " byes" << endl;
206 bool BamStandardIndex::CheckMagicNumber(void) {
208 // check 'magic number' to see if file is BAI index
210 size_t elementsRead = fread(magic, sizeof(char), 4, m_indexStream);
211 if ( elementsRead != 4 ) {
212 cerr << "BamStandardIndex ERROR: could not read format 'magic number'" << endl;
216 // compare to expected value
217 if ( strncmp(magic, BamStandardIndex::BAI_MAGIC, 4) != 0 ) {
218 cerr << "BamStandardIndex ERROR: invalid format" << endl;
226 void BamStandardIndex::ClearReferenceEntry(BaiReferenceEntry& refEntry) {
228 refEntry.Bins.clear();
229 refEntry.LinearOffsets.clear();
232 void BamStandardIndex::CloseFile(void) {
236 fclose(m_indexStream);
238 // clear index file summary data
239 m_indexFileSummary.clear();
241 // clean up I/O buffer
247 // builds index from associated BAM file & writes out to index file
248 bool BamStandardIndex::Create(void) {
250 // return false if BamReader is invalid or not open
251 if ( m_reader == 0 || !m_reader->IsOpen() ) {
252 cerr << "BamStandardIndex ERROR: BamReader is not open"
253 << ", aborting index creation" << endl;
258 if ( !m_reader->Rewind() ) {
259 cerr << "BamStandardIndex ERROR: could not rewind BamReader to create index"
260 << ", aborting index creation" << endl;
264 // open new index file (read & write)
265 string indexFilename = m_reader->Filename() + Extension();
266 if ( !OpenFile(indexFilename, "w+b") ) {
267 cerr << "BamStandardIndex ERROR: could not open ouput index file: " << indexFilename
268 << ", aborting index creation" << endl;
272 // initialize BaiFileSummary with number of references
273 const int& numReferences = m_reader->GetReferenceCount();
274 ReserveForSummary(numReferences);
276 // initialize output file
277 bool createdOk = true;
278 createdOk &= WriteHeader();
280 // set up bin, ID, offset, & coordinate markers
281 const uint32_t defaultValue = 0xffffffffu;
282 uint32_t currentBin = defaultValue;
283 uint32_t lastBin = defaultValue;
284 int32_t currentRefID = defaultValue;
285 int32_t lastRefID = defaultValue;
286 uint64_t currentOffset = (uint64_t)m_reader->Tell();
287 uint64_t lastOffset = currentOffset;
288 int32_t lastPosition = defaultValue;
290 // iterate through alignments in BAM file
292 BaiReferenceEntry refEntry;
293 while ( m_reader->LoadNextAlignment(al) ) {
295 // changed to new reference
296 if ( lastRefID != al.RefID ) {
298 // if not first reference, save previous reference data
299 if ( lastRefID != (int32_t)defaultValue ) {
301 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
302 createdOk &= WriteReferenceEntry(refEntry);
303 ClearReferenceEntry(refEntry);
305 // write any empty references between (but *NOT* including) lastRefID & al.RefID
306 for ( int i = lastRefID+1; i < al.RefID; ++i ) {
307 BaiReferenceEntry emptyEntry(i);
308 createdOk &= WriteReferenceEntry(emptyEntry);
311 // update bin markers
312 currentOffset = lastOffset;
315 currentRefID = al.RefID;
319 // write any empty references up to (but *NOT* including) al.RefID
321 for ( int i = 0; i < al.RefID; ++i ) {
322 BaiReferenceEntry emptyEntry(i);
323 createdOk &= WriteReferenceEntry(emptyEntry);
327 // update reference markers
328 refEntry.ID = al.RefID;
329 lastRefID = al.RefID;
330 lastBin = defaultValue;
333 // if lastPosition greater than current alignment position - file not sorted properly
334 else if ( lastPosition > al.Position ) {
335 cerr << "BamStandardIndex ERROR: BAM file is not properly sorted by coordinate"
336 << ", aborting index creation"
338 << "At alignment: " << al.Name
339 << " : previous position " << lastPosition
340 << " > this alignment position " << al.Position
341 << " on reference id: " << al.RefID << endl;
345 // if alignment's ref ID is valid & its bin is not a 'leaf'
346 if ( (al.RefID >= 0) && (al.Bin < 4681) )
347 SaveLinearOffsetEntry(refEntry.LinearOffsets, al.Position, al.GetEndPosition(), lastOffset);
349 // changed to new BAI bin
350 if ( al.Bin != lastBin ) {
352 // if not first bin on reference, save previous bin data
353 if ( currentBin != defaultValue )
354 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
357 currentOffset = lastOffset;
360 currentRefID = al.RefID;
362 // if invalid RefID, break out
363 if ( currentRefID < 0 )
367 // make sure that current file pointer is beyond lastOffset
368 if ( m_reader->Tell() <= (int64_t)lastOffset ) {
369 cerr << "BamStandardIndex ERROR: calculating offsets failed"
370 << ", aborting index creation" << endl;
374 // update lastOffset & lastPosition
375 lastOffset = m_reader->Tell();
376 lastPosition = al.Position;
379 // after finishing alignments, if any data was read, check:
380 if ( currentRefID >= 0 ) {
382 // store last alignment chunk to its bin, then write last reference entry with data
383 SaveAlignmentChunkToBin(refEntry.Bins, currentBin, currentOffset, lastOffset);
384 createdOk &= WriteReferenceEntry(refEntry);
386 // then write any empty references remaining at end of file
387 for ( int i = currentRefID+1; i < numReferences; ++i ) {
388 BaiReferenceEntry emptyEntry(i);
389 createdOk &= WriteReferenceEntry(emptyEntry);
393 // rewind reader now that we're done building
394 createdOk &= m_reader->Rewind();
400 // returns format's file extension
401 const string BamStandardIndex::Extension(void) {
402 return BamStandardIndex::BAI_EXTENSION;
405 bool BamStandardIndex::GetOffsets(const BamRegion& region, vector<int64_t>& offsets) {
407 // cannot calculate offsets if unknown/invalid reference ID requested
408 if ( region.LeftRefID < 0 || region.LeftRefID >= (int)m_indexFileSummary.size() )
411 // retrieve index summary for left bound reference
412 const BaiReferenceSummary& refSummary = m_indexFileSummary.at(region.LeftRefID);
414 // set up region boundaries based on actual BamReader data
417 if ( !AdjustRegion(region, begin, end) )
420 // retrieve all candidate bin IDs for region
421 set<uint16_t> candidateBins;
422 CalculateCandidateBins(begin, end, candidateBins);
424 // use reference's linear offsets to calculate the minimum offset
425 // that must be considered to find overlap
426 const uint64_t& minOffset = CalculateMinOffset(refSummary, begin);
428 // attempt to use reference summary, minOffset, & candidateBins to calculate offsets
429 if ( !CalculateCandidateOffsets(refSummary, minOffset, candidateBins, offsets) )
432 // ensure that offsets are sorted before returning
433 sort( offsets.begin(), offsets.end() );
439 // returns whether reference has alignments or no
440 bool BamStandardIndex::HasAlignments(const int& referenceID) const {
441 if ( referenceID < 0 || referenceID >= (int)m_indexFileSummary.size() )
443 const BaiReferenceSummary& refSummary = m_indexFileSummary.at(referenceID);
444 return ( refSummary.NumBins > 0 );
447 bool BamStandardIndex::IsFileOpen(void) const {
448 return ( m_indexStream != 0 );
451 // attempts to use index data to jump to @region, returns success/fail
452 // a "successful" jump indicates no error, but not whether this region has data
453 // * thus, the method sets a flag to indicate whether there are alignments
454 // available after the jump position
455 bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
457 // skip if reader is not valid or is not open
458 if ( m_reader == 0 || !m_reader->IsOpen() )
461 // calculate offsets for this region
462 // if failed, print message, set flag, and return failure
463 vector<int64_t> offsets;
464 if ( !GetOffsets(region, offsets) ) {
465 cerr << "BamStandardIndex ERROR: could not jump"
466 << ", unable to retrieve offsets for region" << endl;
467 *hasAlignmentsInRegion = false;
471 // if no offsets retrieved, set flag
472 if ( offsets.empty() )
473 *hasAlignmentsInRegion = false;
475 // iterate through candidate offsets
478 vector<int64_t>::const_iterator offsetIter = offsets.begin();
479 vector<int64_t>::const_iterator offsetEnd = offsets.end();
480 for ( ; offsetIter != offsetEnd; ++offsetIter) {
482 // attempt seek & load first available alignment
483 // set flag to true if data exists
484 result &= m_reader->Seek(*offsetIter);
485 *hasAlignmentsInRegion = m_reader->LoadNextAlignment(al);
487 // if this alignment corresponds to desired position
488 // return success of seeking back to the offset before the 'current offset' (to cover overlaps)
489 if ( ((al.RefID == region.LeftRefID) &&
490 ((al.Position + al.Length) > region.LeftPosition)) ||
491 (al.RefID > region.LeftRefID) )
493 if ( offsetIter != offsets.begin() )
495 return m_reader->Seek(*offsetIter);
499 // if error in jumping, print message & set flag
501 cerr << "BamStandardIndex ERROR: could not jump"
502 << ", there was a problem seeking in BAM file" << endl;
503 *hasAlignmentsInRegion = false;
506 // return success/failure
510 // loads existing data from file into memory
511 bool BamStandardIndex::Load(const std::string& filename) {
513 // attempt open index file (read-only)
514 if ( !OpenFile(filename, "rb") ) {
515 cerr << "BamStandardIndex ERROR: could not open input index file: " << filename
516 << ", aborting index load" << endl;
520 // if invalid format 'magic number', close & return failure
521 if ( !CheckMagicNumber() ) {
522 cerr << "BamStandardIndex ERROR: unexpected format for index file: " << filename
523 << ", aborting index load" << endl;
528 // attempt to load index file summary, return success/failure
529 if ( !SummarizeIndexFile() ) {
530 cerr << "BamStandardIndex ERROR: could not generate a summary of index file " << filename
531 << ", aborting index load" << endl;
536 // if we get here, index summary is loaded OK
540 uint64_t BamStandardIndex::LookupLinearOffset(const BaiReferenceSummary& refSummary, const int& index) {
542 // attempt seek to proper index file position
543 const int64_t linearOffsetFilePosition = (int64_t)refSummary.FirstLinearOffsetFilePosition +
544 index*BamStandardIndex::SIZEOF_LINEAROFFSET;
545 if ( !Seek(linearOffsetFilePosition, SEEK_SET) )
548 // read linear offset from BAI file
549 uint64_t linearOffset(0);
550 if ( !ReadLinearOffset(linearOffset) )
555 void BamStandardIndex::MergeAlignmentChunks(BaiAlignmentChunkVector& chunks) {
557 // skip if chunks are empty, nothing to merge
558 if ( chunks.empty() )
561 // set up merged alignment chunk container
562 BaiAlignmentChunkVector mergedChunks;
563 mergedChunks.push_back( chunks[0] );
565 // iterate over chunks
567 BaiAlignmentChunkVector::iterator chunkIter = chunks.begin();
568 BaiAlignmentChunkVector::iterator chunkEnd = chunks.end();
569 for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {
571 // get 'currentMergeChunk' based on numeric index
572 BaiAlignmentChunk& currentMergeChunk = mergedChunks[i];
574 // get sourceChunk based on source vector iterator
575 BaiAlignmentChunk& sourceChunk = (*chunkIter);
577 // if currentMergeChunk ends where sourceChunk starts, then merge the two
578 if ( currentMergeChunk.Stop>>16 == sourceChunk.Start>>16 )
579 currentMergeChunk.Stop = sourceChunk.Stop;
583 // append sourceChunk after currentMergeChunk
584 mergedChunks.push_back(sourceChunk);
586 // update i, so the next iteration will consider the
587 // recently-appended sourceChunk as new mergeChunk candidate
592 // saved newly-merged chunks into (parameter) chunks
593 chunks = mergedChunks;
596 bool BamStandardIndex::OpenFile(const std::string& filename, const char* mode) {
598 // make sure any previous index file is closed
601 // attempt to open file
602 m_indexStream = fopen(filename.c_str(), mode);
606 bool BamStandardIndex::ReadBinID(uint32_t& binId) {
607 size_t elementsRead = 0;
608 elementsRead += fread(&binId, sizeof(binId), 1, m_indexStream);
609 if ( m_isBigEndian ) SwapEndian_32(binId);
610 return ( elementsRead == 1 );
613 bool BamStandardIndex::ReadBinIntoBuffer(uint32_t& binId, int32_t& numAlignmentChunks) {
618 readOk &= ReadBinID(binId);
619 readOk &= ReadNumAlignmentChunks(numAlignmentChunks);
622 const unsigned int bytesRequested = numAlignmentChunks*BamStandardIndex::SIZEOF_ALIGNMENTCHUNK;
623 readOk &= ReadIntoBuffer(bytesRequested);
625 // return success/failure
629 bool BamStandardIndex::ReadIntoBuffer(const unsigned int& bytesRequested) {
631 // ensure that our buffer is big enough for request
632 BamStandardIndex::CheckBufferSize(m_buffer, m_bufferLength, bytesRequested);
634 // read from BAI file stream
635 size_t bytesRead = fread( m_buffer, sizeof(char), bytesRequested, m_indexStream );
636 return ( bytesRead == (size_t)bytesRequested );
639 bool BamStandardIndex::ReadLinearOffset(uint64_t& linearOffset) {
640 size_t elementsRead = 0;
641 elementsRead += fread(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
642 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
643 return ( elementsRead == 1 );
646 bool BamStandardIndex::ReadNumAlignmentChunks(int& numAlignmentChunks) {
647 size_t elementsRead = 0;
648 elementsRead += fread(&numAlignmentChunks, sizeof(numAlignmentChunks), 1, m_indexStream);
649 if ( m_isBigEndian ) SwapEndian_32(numAlignmentChunks);
650 return ( elementsRead == 1 );
653 bool BamStandardIndex::ReadNumBins(int& numBins) {
654 size_t elementsRead = 0;
655 elementsRead += fread(&numBins, sizeof(numBins), 1, m_indexStream);
656 if ( m_isBigEndian ) SwapEndian_32(numBins);
657 return ( elementsRead == 1 );
660 bool BamStandardIndex::ReadNumLinearOffsets(int& numLinearOffsets) {
661 size_t elementsRead = 0;
662 elementsRead += fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, m_indexStream);
663 if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets);
664 return ( elementsRead == 1 );
667 bool BamStandardIndex::ReadNumReferences(int& numReferences) {
668 size_t elementsRead = 0;
669 elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
670 if ( m_isBigEndian ) SwapEndian_32(numReferences);
671 return ( elementsRead == 1 );
674 void BamStandardIndex::ReserveForSummary(const int& numReferences) {
675 m_indexFileSummary.clear();
676 m_indexFileSummary.assign( numReferences, BaiReferenceSummary() );
679 void BamStandardIndex::SaveAlignmentChunkToBin(BaiBinMap& binMap,
680 const uint32_t& currentBin,
681 const uint64_t& currentOffset,
682 const uint64_t& lastOffset)
684 // create new alignment chunk
685 BaiAlignmentChunk newChunk(currentOffset, lastOffset);
687 // if no entry exists yet for this bin, create one and store alignment chunk
688 BaiBinMap::iterator binIter = binMap.find(currentBin);
689 if ( binIter == binMap.end() ) {
690 BaiAlignmentChunkVector newChunks;
691 newChunks.push_back(newChunk);
692 binMap.insert( pair<uint32_t, BaiAlignmentChunkVector>(currentBin, newChunks));
695 // otherwise, just append alignment chunk
697 BaiAlignmentChunkVector& binChunks = (*binIter).second;
698 binChunks.push_back( newChunk );
702 void BamStandardIndex::SaveBinsSummary(const int& refId, const int& numBins) {
703 BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
704 refSummary.NumBins = numBins;
705 refSummary.FirstBinFilePosition = Tell();
708 void BamStandardIndex::SaveLinearOffsetEntry(BaiLinearOffsetVector& offsets,
709 const int& alignmentStartPosition,
710 const int& alignmentStopPosition,
711 const uint64_t& lastOffset)
713 // get converted offsets
714 const int beginOffset = alignmentStartPosition >> BamStandardIndex::BAM_LIDX_SHIFT;
715 const int endOffset = (alignmentStopPosition - 1) >> BamStandardIndex::BAM_LIDX_SHIFT;
717 // resize vector if necessary
718 int oldSize = offsets.size();
719 int newSize = endOffset + 1;
720 if ( oldSize < newSize )
721 offsets.resize(newSize, 0);
724 for( int i = beginOffset + 1; i <= endOffset; ++i ) {
725 if ( offsets[i] == 0 )
726 offsets[i] = lastOffset;
730 void BamStandardIndex::SaveLinearOffsetsSummary(const int& refId, const int& numLinearOffsets) {
731 BaiReferenceSummary& refSummary = m_indexFileSummary.at(refId);
732 refSummary.NumLinearOffsets = numLinearOffsets;
733 refSummary.FirstLinearOffsetFilePosition = Tell();
736 // seek to position in index file stream
737 bool BamStandardIndex::Seek(const int64_t& position, const int& origin) {
738 return ( fseek64(m_indexStream, position, origin) == 0 );
741 // change the index caching behavior
742 void BamStandardIndex::SetCacheMode(const BamIndex::IndexCacheMode& mode) {
744 // do nothing else here ? cache mode will be ignored from now on, most likely
747 bool BamStandardIndex::SkipBins(const int& numBins) {
749 int32_t numAlignmentChunks;
750 bool skippedOk = true;
751 for (int i = 0; i < numBins; ++i)
752 skippedOk &= ReadBinIntoBuffer(binId, numAlignmentChunks); // results & buffer ignored
756 bool BamStandardIndex::SkipLinearOffsets(const int& numLinearOffsets) {
757 const unsigned int bytesRequested = numLinearOffsets*BamStandardIndex::SIZEOF_LINEAROFFSET;
758 return ReadIntoBuffer(bytesRequested);
761 void BamStandardIndex::SortLinearOffsets(BaiLinearOffsetVector& linearOffsets) {
762 sort( linearOffsets.begin(), linearOffsets.end() );
765 bool BamStandardIndex::SummarizeBins(BaiReferenceSummary& refSummary) {
767 // load number of bins
769 if ( !ReadNumBins(numBins) )
772 // store bins summary for this reference
773 refSummary.NumBins = numBins;
774 refSummary.FirstBinFilePosition = Tell();
776 // attempt skip reference bins, return success/failure
777 if ( !SkipBins(numBins) )
780 // if we get here, bin summarized OK
784 bool BamStandardIndex::SummarizeIndexFile(void) {
786 // load number of reference sequences
788 if ( !ReadNumReferences(numReferences) )
791 // initialize file summary data
792 ReserveForSummary(numReferences);
794 // iterate over reference entries
795 bool loadedOk = true;
796 BaiFileSummary::iterator summaryIter = m_indexFileSummary.begin();
797 BaiFileSummary::iterator summaryEnd = m_indexFileSummary.end();
798 for ( int i = 0; summaryIter != summaryEnd; ++summaryIter, ++i )
799 loadedOk &= SummarizeReference(*summaryIter);
805 bool BamStandardIndex::SummarizeLinearOffsets(BaiReferenceSummary& refSummary) {
807 // load number of linear offsets
808 int numLinearOffsets;
809 if ( !ReadNumLinearOffsets(numLinearOffsets) )
812 // store bin summary data for this reference
813 refSummary.NumLinearOffsets = numLinearOffsets;
814 refSummary.FirstLinearOffsetFilePosition = Tell();
816 // skip linear offsets in index file
817 if ( !SkipLinearOffsets(numLinearOffsets) )
820 // if get here, linear offsets summarized OK
824 bool BamStandardIndex::SummarizeReference(BaiReferenceSummary& refSummary) {
825 bool loadedOk = true;
826 loadedOk &= SummarizeBins(refSummary);
827 loadedOk &= SummarizeLinearOffsets(refSummary);
831 // return position of file pointer in index file stream
832 int64_t BamStandardIndex::Tell(void) const {
833 return ftell64(m_indexStream);
836 bool BamStandardIndex::WriteAlignmentChunk(const BaiAlignmentChunk& chunk) {
838 size_t elementsWritten = 0;
840 // localize alignment chunk offsets
841 uint64_t start = chunk.Start;
842 uint64_t stop = chunk.Stop;
844 // swap endian-ness if necessary
845 if ( m_isBigEndian ) {
846 SwapEndian_64(start);
850 // write to index file
851 elementsWritten += fwrite(&start, sizeof(start), 1, m_indexStream);
852 elementsWritten += fwrite(&stop, sizeof(stop), 1, m_indexStream);
854 // return success/failure of write
855 return ( elementsWritten == 2 );
858 bool BamStandardIndex::WriteAlignmentChunks(BaiAlignmentChunkVector& chunks) {
860 // make sure chunks are merged (simplified) before writing & saving summary
861 MergeAlignmentChunks(chunks);
863 size_t elementsWritten = 0;
866 int32_t chunkCount = chunks.size();
867 if ( m_isBigEndian ) SwapEndian_32(chunkCount);
868 elementsWritten += fwrite(&chunkCount, sizeof(chunkCount), 1, m_indexStream);
870 // iterate over chunks
871 bool chunksOk = true;
872 BaiAlignmentChunkVector::const_iterator chunkIter = chunks.begin();
873 BaiAlignmentChunkVector::const_iterator chunkEnd = chunks.end();
874 for ( ; chunkIter != chunkEnd; ++chunkIter )
875 chunksOk &= WriteAlignmentChunk( (*chunkIter) );
877 // return success/failure of write
878 return ( (elementsWritten == 1) && chunksOk );
881 bool BamStandardIndex::WriteBin(const uint32_t& binId, BaiAlignmentChunkVector& chunks) {
883 size_t elementsWritten = 0;
886 uint32_t binKey = binId;
887 if ( m_isBigEndian ) SwapEndian_32(binKey);
888 elementsWritten += fwrite(&binKey, sizeof(binKey), 1, m_indexStream);
890 // write bin's alignment chunks
891 bool chunksOk = WriteAlignmentChunks(chunks);
893 // return success/failure of write
894 return ( (elementsWritten == 1) && chunksOk );
897 bool BamStandardIndex::WriteBins(const int& refId, BaiBinMap& bins) {
899 size_t elementsWritten = 0;
901 // write number of bins
902 int32_t binCount = bins.size();
903 if ( m_isBigEndian ) SwapEndian_32(binCount);
904 elementsWritten += fwrite(&binCount, sizeof(binCount), 1, m_indexStream);
906 // save summary for reference's bins
907 SaveBinsSummary(refId, bins.size());
911 BaiBinMap::iterator binIter = bins.begin();
912 BaiBinMap::iterator binEnd = bins.end();
913 for ( ; binIter != binEnd; ++binIter )
914 binsOk &= WriteBin( (*binIter).first, (*binIter).second );
916 // return success/failure of write
917 return ( (elementsWritten == 1) && binsOk );
920 bool BamStandardIndex::WriteHeader(void) {
922 size_t elementsWritten = 0;
924 // write magic number
925 elementsWritten += fwrite(BamStandardIndex::BAI_MAGIC, sizeof(char), 4, m_indexStream);
927 // write number of reference sequences
928 int32_t numReferences = m_indexFileSummary.size();
929 if ( m_isBigEndian ) SwapEndian_32(numReferences);
930 elementsWritten += fwrite(&numReferences, sizeof(numReferences), 1, m_indexStream);
932 // return success/failure of write
933 return (elementsWritten == 5);
936 bool BamStandardIndex::WriteLinearOffsets(const int& refId, BaiLinearOffsetVector& linearOffsets) {
938 // make sure linear offsets are sorted before writing & saving summary
939 SortLinearOffsets(linearOffsets);
941 size_t elementsWritten = 0;
943 // write number of linear offsets
944 int32_t offsetCount = linearOffsets.size();
945 if ( m_isBigEndian ) SwapEndian_32(offsetCount);
946 elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, m_indexStream);
948 // save summary for reference's linear offsets
949 SaveLinearOffsetsSummary(refId, linearOffsets.size());
951 // iterate over linear offsets
952 BaiLinearOffsetVector::const_iterator offsetIter = linearOffsets.begin();
953 BaiLinearOffsetVector::const_iterator offsetEnd = linearOffsets.end();
954 for ( ; offsetIter != offsetEnd; ++offsetIter ) {
956 // write linear offset
957 uint64_t linearOffset = (*offsetIter);
958 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
959 elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
962 // return success/failure of write
963 return ( elementsWritten == (size_t)(linearOffsets.size() + 1) );
966 bool BamStandardIndex::WriteReferenceEntry(BaiReferenceEntry& refEntry) {
968 refOk &= WriteBins(refEntry.ID, refEntry.Bins);
969 refOk &= WriteLinearOffsets(refEntry.ID, refEntry.LinearOffsets);