1 // ***************************************************************************
2 // BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg
3 // Marth Lab, Department of Biology, Boston College
4 // All rights reserved.
5 // ---------------------------------------------------------------------------
6 // Last modified: 13 October 2010 (DB)
7 // ---------------------------------------------------------------------------
8 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad
10 // ---------------------------------------------------------------------------
11 // Provides the basic functionality for reading BAM files
12 // ***************************************************************************
21 #include "BamReader.h"
22 using namespace BamTools;
25 struct BamReader::BamReaderPrivate {
27 // -------------------------------
28 // structs, enums, typedefs
29 // -------------------------------
30 enum RegionState { BEFORE_REGION = 0
35 // -------------------------------
37 // -------------------------------
45 int64_t AlignmentsBeginOffset;
50 BamIndex::BamIndexCacheMode IndexCacheMode;
55 // user-specified region values
57 bool HasAlignmentsInRegion;
62 // BAM character constants
63 const char* DNA_LOOKUP;
64 const char* CIGAR_LOOKUP;
66 // constructor & destructor
67 BamReaderPrivate(BamReader* parent);
68 ~BamReaderPrivate(void);
70 // -------------------------------
75 bool Open(const std::string& filename,
76 const std::string& indexFilename,
77 const bool lookForIndex,
78 const bool preferStandardIndex);
80 bool SetRegion(const BamRegion& region);
82 // access alignment data
83 bool GetNextAlignment(BamAlignment& bAlignment);
84 bool GetNextAlignmentCore(BamAlignment& bAlignment);
86 // access auxiliary data
87 int GetReferenceID(const string& refName) const;
90 bool CreateIndex(bool useStandardIndex);
91 void SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode);
97 // ---------------------------------------
98 // reading alignments and auxiliary data
100 // adjusts requested region if necessary (depending on where data actually begins)
101 void AdjustRegion(BamRegion& region);
102 // fills out character data for BamAlignment data
103 bool BuildCharData(BamAlignment& bAlignment);
104 // checks to see if alignment overlaps current region
105 RegionState IsOverlap(BamAlignment& bAlignment);
106 // retrieves header text from BAM file
107 void LoadHeaderData(void);
108 // retrieves BAM alignment under file pointer
109 bool LoadNextAlignment(BamAlignment& bAlignment);
110 // builds reference data structure from BAM file
111 void LoadReferenceData(void);
112 // mark references with 'HasAlignments' status
113 void MarkReferences(void);
115 // ---------------------------------
116 // index file handling
118 // clear out inernal index data structure
119 void ClearIndex(void);
120 // loads index from BAM index file
121 bool LoadIndex(const bool lookForIndex, const bool preferStandardIndex);
124 // -----------------------------------------------------
125 // BamReader implementation (wrapper around BRPrivate)
126 // -----------------------------------------------------
128 BamReader::BamReader(void) {
129 d = new BamReaderPrivate(this);
133 BamReader::~BamReader(void) {
139 void BamReader::Close(void) { d->Close(); }
140 bool BamReader::HasIndex(void) const { return d->HasIndex; }
141 bool BamReader::IsIndexLoaded(void) const { return HasIndex(); }
142 bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }
143 bool BamReader::Jump(int refID, int position) { return d->SetRegion( BamRegion(refID, position) ); }
144 bool BamReader::Open(const std::string& filename,
145 const std::string& indexFilename,
146 const bool lookForIndex,
147 const bool preferStandardIndex)
149 return d->Open(filename, indexFilename, lookForIndex, preferStandardIndex);
151 bool BamReader::Rewind(void) { return d->Rewind(); }
152 bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }
153 bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {
154 return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );
157 // access alignment data
158 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }
159 bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }
161 // access auxiliary data
162 const string BamReader::GetHeaderText(void) const { return d->HeaderText; }
163 int BamReader::GetReferenceCount(void) const { return d->References.size(); }
164 const RefVector& BamReader::GetReferenceData(void) const { return d->References; }
165 int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }
166 const std::string BamReader::GetFilename(void) const { return d->Filename; }
169 bool BamReader::CreateIndex(bool useStandardIndex) { return d->CreateIndex(useStandardIndex); }
170 void BamReader::SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode) { d->SetIndexCacheMode(mode); }
172 // -----------------------------------------------------
173 // BamReaderPrivate implementation
174 // -----------------------------------------------------
177 BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)
180 , AlignmentsBeginOffset(0)
181 , IndexCacheMode(BamIndex::LimitedIndexCaching)
182 , HasAlignmentsInRegion(true)
184 , DNA_LOOKUP("=ACMGRSVTWYHKDBN")
185 , CIGAR_LOOKUP("MIDNSHP")
187 IsBigEndian = SystemIsBigEndian();
191 BamReader::BamReaderPrivate::~BamReaderPrivate(void) {
195 // adjusts requested region if necessary (depending on where data actually begins)
196 void BamReader::BamReaderPrivate::AdjustRegion(BamRegion& region) {
198 // check for valid index first
199 if ( Index == 0 ) return;
201 // see if any references in region have alignments
202 HasAlignmentsInRegion = false;
203 int currentId = region.LeftRefID;
204 while ( currentId <= region.RightRefID ) {
205 HasAlignmentsInRegion = Index->HasAlignments(currentId);
206 if ( HasAlignmentsInRegion ) break;
210 // if no data found on any reference in region
211 if ( !HasAlignmentsInRegion ) return;
213 // if left bound of desired region had no data, use first reference that had data
214 // otherwise, leave requested region as-is
215 if ( currentId != region.LeftRefID ) {
216 region.LeftRefID = currentId;
217 region.LeftPosition = 0;
221 // fills out character data for BamAlignment data
222 bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {
224 // calculate character lengths/offsets
225 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
226 const unsigned int seqDataOffset = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);
227 const unsigned int qualDataOffset = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;
228 const unsigned int tagDataOffset = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;
229 const unsigned int tagDataLength = dataLength - tagDataOffset;
231 // set up char buffers
232 const char* allCharData = bAlignment.SupportData.AllCharData.data();
233 const char* seqData = ((const char*)allCharData) + seqDataOffset;
234 const char* qualData = ((const char*)allCharData) + qualDataOffset;
235 char* tagData = ((char*)allCharData) + tagDataOffset;
237 // store alignment name (relies on null char in name as terminator)
238 bAlignment.Name.assign((const char*)(allCharData));
240 // save query sequence
241 bAlignment.QueryBases.clear();
242 bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);
243 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
244 char singleBase = DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
245 bAlignment.QueryBases.append(1, singleBase);
248 // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
249 bAlignment.Qualities.clear();
250 bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);
251 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
252 char singleQuality = (char)(qualData[i]+33);
253 bAlignment.Qualities.append(1, singleQuality);
256 // if QueryBases is empty (and this is a allowed case)
257 if ( bAlignment.QueryBases.empty() )
258 bAlignment.AlignedBases = bAlignment.QueryBases;
260 // if QueryBases contains data, then build AlignedBases using CIGAR data
263 // resize AlignedBases
264 bAlignment.AlignedBases.clear();
265 bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);
267 // iterate over CigarOps
269 vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();
270 vector<CigarOp>::const_iterator cigarEnd = bAlignment.CigarData.end();
271 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
273 const CigarOp& op = (*cigarIter);
278 bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases
282 k += op.Length; // for 'S' - soft clip, skip over query bases
286 bAlignment.AlignedBases.append(op.Length, '-'); // for 'D' - write gap character
290 bAlignment.AlignedBases.append( op.Length, '*' ); // for 'P' - write padding character
294 bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in original query sequence
298 break; // for 'H' - hard clip, do nothing to AlignedBases, move to next op
301 fprintf(stderr, "ERROR: Invalid Cigar op type\n"); // shouldn't get here
307 // -----------------------
308 // Added: 3-25-2010 DB
309 // Fixed: endian-correctness for tag data
310 // -----------------------
313 while ( (unsigned int)i < tagDataLength ) {
315 i += 2; // skip tag type (e.g. "RG", "NM", etc)
316 uint8_t type = toupper(tagData[i]); // lower & upper case letters have same meaning
317 ++i; // skip value type
327 SwapEndian_16p(&tagData[i]);
328 i += sizeof(uint16_t);
333 SwapEndian_32p(&tagData[i]);
334 i += sizeof(uint32_t);
338 SwapEndian_64p(&tagData[i]);
339 i += sizeof(uint64_t);
344 while (tagData[i]) { ++i; }
345 ++i; // increment one more for null terminator
349 fprintf(stderr, "ERROR: Invalid tag value type\n"); // shouldn't get here
356 bAlignment.TagData.clear();
357 bAlignment.TagData.resize(tagDataLength);
358 memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);
360 // clear the core-only flag
361 bAlignment.SupportData.HasCoreOnly = false;
367 // clear index data structure
368 void BamReader::BamReaderPrivate::ClearIndex(void) {
374 // closes the BAM file
375 void BamReader::BamReaderPrivate::Close(void) {
377 // close BGZF file stream
380 // clear out index data
383 // clear out header data
386 // clear out region flags
390 // creates index for BAM file, saves to file
391 // default behavior is to create the BAM standard index (".bai")
392 // set flag to false to create the BamTools-specific index (".bti")
393 bool BamReader::BamReaderPrivate::CreateIndex(bool useStandardIndex) {
395 // clear out prior index data
398 // create index based on type requested
399 if ( useStandardIndex )
400 Index = new BamStandardIndex(&mBGZF, Parent);
402 Index = new BamToolsIndex(&mBGZF, Parent);
404 // set index cache mode to full for writing
405 Index->SetCacheMode(BamIndex::FullIndexCaching);
409 ok &= Index->Build();
412 // mark empty references
415 // attempt to save index data to file
416 ok &= Index->Write(Filename);
418 // set client's desired index cache mode
419 Index->SetCacheMode(IndexCacheMode);
421 // return success/fail of both building & writing index
425 // get next alignment (from specified region, if given)
426 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {
428 // if valid alignment found, attempt to parse char data, and return success/failure
429 if ( GetNextAlignmentCore(bAlignment) )
430 return BuildCharData(bAlignment);
432 // no valid alignment found
436 // retrieves next available alignment core data (returns success/fail)
437 // ** DOES NOT parse any character data (read name, bases, qualities, tag data)
438 // these can be accessed, if necessary, from the supportData
439 // useful for operations requiring ONLY positional or other alignment-related information
440 bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {
442 // if region is set but has no alignments
443 if ( !Region.isNull() && !HasAlignmentsInRegion )
446 // if valid alignment available
447 if ( LoadNextAlignment(bAlignment) ) {
449 // set core-only flag
450 bAlignment.SupportData.HasCoreOnly = true;
452 // if region not specified with at least a left boundary, return success
453 if ( !Region.isLeftBoundSpecified() ) return true;
455 // determine region state (before, within, after)
456 BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);
458 // if alignment lies after region, return false
459 if ( state == AFTER_REGION ) return false;
461 while ( state != WITHIN_REGION ) {
462 // if no valid alignment available (likely EOF) return failure
463 if ( !LoadNextAlignment(bAlignment) ) return false;
464 // if alignment lies after region, return false (no available read within region)
465 state = IsOverlap(bAlignment);
466 if ( state == AFTER_REGION ) return false;
469 // return success (alignment found that overlaps region)
473 // no valid alignment
477 // returns RefID for given RefName (returns References.size() if not found)
478 int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {
480 // retrieve names from reference data
481 vector<string> refNames;
482 RefVector::const_iterator refIter = References.begin();
483 RefVector::const_iterator refEnd = References.end();
484 for ( ; refIter != refEnd; ++refIter)
485 refNames.push_back( (*refIter).RefName );
487 // return 'index-of' refName ( if not found, returns refNames.size() )
488 return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
491 // returns region state - whether alignment ends before, overlaps, or starts after currently specified region
492 // this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
493 BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {
495 // if alignment is on any reference sequence before left bound
496 if ( bAlignment.RefID < Region.LeftRefID ) return BEFORE_REGION;
498 // if alignment starts on left bound reference
499 else if ( bAlignment.RefID == Region.LeftRefID ) {
501 // if alignment starts at or after left boundary
502 if ( bAlignment.Position >= Region.LeftPosition) {
504 // if right boundary is specified AND
505 // left/right boundaries are on same reference AND
506 // alignment starts past right boundary
507 if ( Region.isRightBoundSpecified() &&
508 Region.LeftRefID == Region.RightRefID &&
509 bAlignment.Position > Region.RightPosition )
512 // otherwise, alignment is within region
513 return WITHIN_REGION;
516 // alignment starts before left boundary
518 // check if alignment overlaps left boundary
519 if ( bAlignment.GetEndPosition() >= Region.LeftPosition ) return WITHIN_REGION;
520 else return BEFORE_REGION;
524 // alignment starts on a reference after the left bound
527 // if region has a right boundary
528 if ( Region.isRightBoundSpecified() ) {
530 // alignment is on reference between boundaries
531 if ( bAlignment.RefID < Region.RightRefID ) return WITHIN_REGION;
533 // alignment is on reference after right boundary
534 else if ( bAlignment.RefID > Region.RightRefID ) return AFTER_REGION;
536 // alignment is on right bound reference
538 // check if alignment starts before or at right boundary
539 if ( bAlignment.Position <= Region.RightPosition ) return WITHIN_REGION;
540 else return AFTER_REGION;
544 // otherwise, alignment is after left bound reference, but there is no right boundary
545 else return WITHIN_REGION;
549 // load BAM header data
550 void BamReader::BamReaderPrivate::LoadHeaderData(void) {
552 // check to see if proper BAM header
554 if (mBGZF.Read(buffer, 4) != 4) {
555 fprintf(stderr, "Could not read header type\n");
559 if (strncmp(buffer, "BAM\001", 4)) {
560 fprintf(stderr, "wrong header type!\n");
564 // get BAM header text length
565 mBGZF.Read(buffer, 4);
566 unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);
567 if ( IsBigEndian ) SwapEndian_32(headerTextLength);
569 // get BAM header text
570 char* headerText = (char*)calloc(headerTextLength + 1, 1);
571 mBGZF.Read(headerText, headerTextLength);
572 HeaderText = (string)((const char*)headerText);
574 // clean up calloc-ed temp variable
578 // load existing index data from BAM index file (".bti" OR ".bai"), return success/fail
579 bool BamReader::BamReaderPrivate::LoadIndex(const bool lookForIndex, const bool preferStandardIndex) {
581 // clear out any existing index data
584 // if no index filename provided, so we need to look for available index files
585 if ( IndexFilename.empty() ) {
587 // attempt to load BamIndex based on current Filename provided & preferStandardIndex flag
588 const BamIndex::PreferredIndexType type = (preferStandardIndex ? BamIndex::STANDARD : BamIndex::BAMTOOLS);
589 Index = BamIndex::FromBamFilename(Filename, &mBGZF, Parent, type);
591 // if null, return failure
592 if ( Index == 0 ) return false;
594 // generate proper IndexFilename based on type of index created
595 IndexFilename = Filename + Index->Extension();
600 // attempt to load BamIndex based on IndexFilename provided by client
601 Index = BamIndex::FromIndexFilename(IndexFilename, &mBGZF, Parent);
603 // if null, return failure
604 if ( Index == 0 ) return false;
607 // set cache mode for BamIndex
608 Index->SetCacheMode(IndexCacheMode);
610 // loading the index data from file
611 HasIndex = Index->Load(IndexFilename);
613 // mark empty references
616 // return index status
620 // populates BamAlignment with alignment data under file pointer, returns success/fail
621 bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {
623 // read in the 'block length' value, make sure it's not zero
625 mBGZF.Read(buffer, 4);
626 bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);
627 if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }
628 if ( bAlignment.SupportData.BlockLength == 0 ) return false;
630 // read in core alignment data, make sure the right size of data was read
631 char x[BAM_CORE_SIZE];
632 if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) return false;
635 for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) )
636 SwapEndian_32p(&x[i]);
639 // set BamAlignment 'core' and 'support' data
640 bAlignment.RefID = BgzfData::UnpackSignedInt(&x[0]);
641 bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);
643 unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);
644 bAlignment.Bin = tempValue >> 16;
645 bAlignment.MapQuality = tempValue >> 8 & 0xff;
646 bAlignment.SupportData.QueryNameLength = tempValue & 0xff;
648 tempValue = BgzfData::UnpackUnsignedInt(&x[12]);
649 bAlignment.AlignmentFlag = tempValue >> 16;
650 bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;
652 bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);
653 bAlignment.MateRefID = BgzfData::UnpackSignedInt(&x[20]);
654 bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);
655 bAlignment.InsertSize = BgzfData::UnpackSignedInt(&x[28]);
657 // set BamAlignment length
658 bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;
660 // read in character data - make sure proper data size was read
661 bool readCharDataOK = false;
662 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
663 char* allCharData = (char*)calloc(sizeof(char), dataLength);
665 if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) {
667 // store 'allCharData' in supportData structure
668 bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
671 readCharDataOK = true;
674 // need to calculate this here so that BamAlignment::GetEndPosition() performs correctly,
675 // even when BamReader::GetNextAlignmentCore() is called
676 const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;
677 uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
679 bAlignment.CigarData.clear();
680 bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);
681 for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {
684 if ( IsBigEndian ) SwapEndian_32(cigarData[i]);
686 // build CigarOp structure
687 op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
688 op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
691 bAlignment.CigarData.push_back(op);
696 return readCharDataOK;
699 // loads reference data from BAM file
700 void BamReader::BamReaderPrivate::LoadReferenceData(void) {
702 // get number of reference sequences
704 mBGZF.Read(buffer, 4);
705 unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);
706 if ( IsBigEndian ) SwapEndian_32(numberRefSeqs);
707 if ( numberRefSeqs == 0 ) return;
708 References.reserve((int)numberRefSeqs);
710 // iterate over all references in header
711 for (unsigned int i = 0; i != numberRefSeqs; ++i) {
713 // get length of reference name
714 mBGZF.Read(buffer, 4);
715 unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);
716 if ( IsBigEndian ) SwapEndian_32(refNameLength);
717 char* refName = (char*)calloc(refNameLength, 1);
719 // get reference name and reference sequence length
720 mBGZF.Read(refName, refNameLength);
721 mBGZF.Read(buffer, 4);
722 int refLength = BgzfData::UnpackSignedInt(buffer);
723 if ( IsBigEndian ) SwapEndian_32(refLength);
725 // store data for reference
727 aReference.RefName = (string)((const char*)refName);
728 aReference.RefLength = refLength;
729 References.push_back(aReference);
731 // clean up calloc-ed temp variable
736 // mark references with no alignment data
737 void BamReader::BamReaderPrivate::MarkReferences(void) {
739 // ensure index is available
740 if ( !HasIndex ) return;
742 // mark empty references
743 for ( int i = 0; i < (int)References.size(); ++i )
744 References.at(i).RefHasAlignments = Index->HasAlignments(i);
747 // opens BAM file (and index)
748 bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename, const bool lookForIndex, const bool preferStandardIndex) {
752 IndexFilename = indexFilename;
754 // open the BGZF file for reading, return false on failure
755 if ( !mBGZF.Open(filename, "rb") ) return false;
757 // retrieve header text & reference data
761 // store file offset of first alignment
762 AlignmentsBeginOffset = mBGZF.Tell();
764 // if no index filename provided
765 if ( IndexFilename.empty() ) {
767 // client did not specify that index SHOULD be found
768 // useful for cases where sequential access is all that is required
769 if ( !lookForIndex ) return true;
771 // otherwise, look for index file, return success/fail
772 return LoadIndex(lookForIndex, preferStandardIndex) ;
775 // client supplied an index filename
776 // attempt to load index data, return success/fail
777 return LoadIndex(lookForIndex, preferStandardIndex);
780 // returns BAM file pointer to beginning of alignment data
781 bool BamReader::BamReaderPrivate::Rewind(void) {
783 // rewind to first alignment, return false if unable to seek
784 if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;
786 // retrieve first alignment data, return false if unable to read
788 if ( !LoadNextAlignment(al) ) return false;
790 // reset default region info using first alignment in file
792 HasAlignmentsInRegion = true;
794 // rewind back to beginning of first alignment
795 // return success/fail of seek
796 return mBGZF.Seek(AlignmentsBeginOffset);
799 // change the index caching behavior
800 void BamReader::BamReaderPrivate::SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode) {
801 IndexCacheMode = mode;
802 if ( Index == 0 ) return;
803 Index->SetCacheMode(mode);
806 // asks Index to attempt a Jump() to specified region
807 // returns success/failure
808 bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {
810 // clear out any prior BamReader region data
812 // N.B. - this is cleared so that BamIndex now has free reign to call
813 // GetNextAlignmentCore() and do overlap checking without worrying about BamReader
814 // performing any overlap checking of its own and moving on to the next read... Calls
815 // to GetNextAlignmentCore() with no Region set, simply return the next alignment.
816 // This ensures that the Index is able to do just that. (All without exposing
817 // LoadNextAlignment() to the public API, and potentially confusing clients with the nomenclature)
820 // check for existing index
821 if ( !HasIndex ) return false;
823 // adjust region if necessary to reflect where data actually begins
824 BamRegion adjustedRegion(region);
825 AdjustRegion(adjustedRegion);
827 // if no data present, return true
828 // not an error, but BamReader knows that no data is there for future alignment access
829 // (this is useful in a MultiBamReader setting where some BAM files may lack data in regions
830 // that other BAMs have data)
831 if ( !HasAlignmentsInRegion ) {
832 Region = adjustedRegion;
836 // attempt jump to user-specified region return false if jump could not be performed at all
837 // (invalid index, unknown reference, etc)
839 // Index::Jump() is allowed to modify the HasAlignmentsInRegion flag
840 // * This covers case where a region is requested that lies beyond the last alignment on a reference
841 // If this occurs, any subsequent calls to GetNexAlignment[Core] simply return false
842 // BamMultiReader is then able to successfully pull alignments from a region from multiple files
843 // even if one or more have no data.
844 if ( !Index->Jump(adjustedRegion, &HasAlignmentsInRegion) ) return false;
846 // save region and return success
847 Region = adjustedRegion;