1 // ***************************************************************************
\r
2 // BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg
\r
3 // Marth Lab, Department of Biology, Boston College
\r
4 // All rights reserved.
\r
5 // ---------------------------------------------------------------------------
\r
6 // Last modified: 3 September 2010 (DB)
\r
7 // ---------------------------------------------------------------------------
\r
8 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad
\r
10 // ---------------------------------------------------------------------------
\r
11 // Provides the basic functionality for reading BAM files
\r
12 // ***************************************************************************
\r
15 #include <algorithm>
\r
21 // BamTools includes
\r
23 #include "BamReader.h"
\r
24 #include "BamIndex.h"
\r
25 using namespace BamTools;
\r
26 using namespace std;
\r
28 struct BamReader::BamReaderPrivate {
\r
30 // -------------------------------
\r
31 // structs, enums, typedefs
\r
32 // -------------------------------
\r
33 enum RegionState { BEFORE_REGION = 0
\r
38 // -------------------------------
\r
40 // -------------------------------
\r
42 // general file data
\r
47 RefVector References;
\r
49 int64_t AlignmentsBeginOffset;
\r
51 string IndexFilename;
\r
56 // user-specified region values
\r
58 bool IsLeftBoundSpecified;
\r
59 bool IsRightBoundSpecified;
\r
61 bool IsRegionSpecified;
\r
68 // BAM character constants
\r
69 const char* DNA_LOOKUP;
\r
70 const char* CIGAR_LOOKUP;
\r
72 // -------------------------------
\r
73 // constructor & destructor
\r
74 // -------------------------------
\r
75 BamReaderPrivate(BamReader* parent);
\r
76 ~BamReaderPrivate(void);
\r
78 // -------------------------------
\r
79 // "public" interface
\r
80 // -------------------------------
\r
84 bool Jump(int refID, int position);
\r
85 bool Open(const std::string& filename,
\r
86 const std::string& indexFilename,
\r
87 const bool lookForIndex,
\r
88 const bool preferStandardIndex);
\r
90 bool SetRegion(const BamRegion& region);
\r
92 // access alignment data
\r
93 bool GetNextAlignment(BamAlignment& bAlignment);
\r
94 bool GetNextAlignmentCore(BamAlignment& bAlignment);
\r
96 // access auxiliary data
\r
97 int GetReferenceID(const string& refName) const;
\r
100 bool CreateIndex(bool useStandardIndex);
\r
102 // -------------------------------
\r
103 // internal methods
\r
104 // -------------------------------
\r
106 // *** reading alignments and auxiliary data *** //
\r
108 // fills out character data for BamAlignment data
\r
109 bool BuildCharData(BamAlignment& bAlignment);
\r
110 // checks to see if alignment overlaps current region
\r
111 RegionState IsOverlap(BamAlignment& bAlignment);
\r
112 // retrieves header text from BAM file
\r
113 void LoadHeaderData(void);
\r
114 // retrieves BAM alignment under file pointer
\r
115 bool LoadNextAlignment(BamAlignment& bAlignment);
\r
116 // builds reference data structure from BAM file
\r
117 void LoadReferenceData(void);
\r
119 // *** index file handling *** //
\r
121 // clear out inernal index data structure
\r
122 void ClearIndex(void);
\r
123 // loads index from BAM index file
\r
124 bool LoadIndex(const bool lookForIndex, const bool preferStandardIndex);
\r
127 // -----------------------------------------------------
\r
128 // BamReader implementation (wrapper around BRPrivate)
\r
129 // -----------------------------------------------------
\r
131 BamReader::BamReader(void) {
\r
132 d = new BamReaderPrivate(this);
\r
136 BamReader::~BamReader(void) {
\r
142 void BamReader::Close(void) { d->Close(); }
\r
143 bool BamReader::IsIndexLoaded(void) const { return d->IsIndexLoaded; }
\r
144 bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }
\r
145 bool BamReader::Jump(int refID, int position)
\r
147 d->Region.LeftRefID = refID;
\r
148 d->Region.LeftPosition = position;
\r
149 d->IsLeftBoundSpecified = true;
\r
150 d->IsRightBoundSpecified = false;
\r
151 return d->Jump(refID, position);
\r
153 bool BamReader::Open(const std::string& filename,
\r
154 const std::string& indexFilename,
\r
155 const bool lookForIndex,
\r
156 const bool preferStandardIndex)
\r
158 return d->Open(filename, indexFilename, lookForIndex, preferStandardIndex);
\r
160 bool BamReader::Rewind(void) { return d->Rewind(); }
\r
161 bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }
\r
162 bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {
\r
163 return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );
\r
166 // access alignment data
\r
167 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }
\r
168 bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }
\r
170 // access auxiliary data
\r
171 const string BamReader::GetHeaderText(void) const { return d->HeaderText; }
\r
172 int BamReader::GetReferenceCount(void) const { return d->References.size(); }
\r
173 const RefVector& BamReader::GetReferenceData(void) const { return d->References; }
\r
174 int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }
\r
175 const std::string BamReader::GetFilename(void) const { return d->Filename; }
\r
177 // index operations
\r
178 bool BamReader::CreateIndex(bool useStandardIndex) { return d->CreateIndex(useStandardIndex); }
\r
180 // -----------------------------------------------------
\r
181 // BamReaderPrivate implementation
\r
182 // -----------------------------------------------------
\r
185 BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)
\r
187 , IsIndexLoaded(false)
\r
188 , AlignmentsBeginOffset(0)
\r
189 , IsLeftBoundSpecified(false)
\r
190 , IsRightBoundSpecified(false)
\r
191 , IsRegionSpecified(false)
\r
195 , DNA_LOOKUP("=ACMGRSVTWYHKDBN")
\r
196 , CIGAR_LOOKUP("MIDNSHP")
\r
198 IsBigEndian = SystemIsBigEndian();
\r
202 BamReader::BamReaderPrivate::~BamReaderPrivate(void) {
\r
206 bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {
\r
208 // calculate character lengths/offsets
\r
209 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
\r
210 const unsigned int seqDataOffset = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);
\r
211 const unsigned int qualDataOffset = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;
\r
212 const unsigned int tagDataOffset = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;
\r
213 const unsigned int tagDataLength = dataLength - tagDataOffset;
\r
215 // set up char buffers
\r
216 const char* allCharData = bAlignment.SupportData.AllCharData.data();
\r
217 const char* seqData = ((const char*)allCharData) + seqDataOffset;
\r
218 const char* qualData = ((const char*)allCharData) + qualDataOffset;
\r
219 char* tagData = ((char*)allCharData) + tagDataOffset;
\r
221 // store alignment name (relies on null char in name as terminator)
\r
222 bAlignment.Name.assign((const char*)(allCharData));
\r
224 // save query sequence
\r
225 bAlignment.QueryBases.clear();
\r
226 bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);
\r
227 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
\r
228 char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
\r
229 bAlignment.QueryBases.append(1, singleBase);
\r
232 // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
\r
233 bAlignment.Qualities.clear();
\r
234 bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);
\r
235 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
\r
236 char singleQuality = (char)(qualData[i]+33);
\r
237 bAlignment.Qualities.append(1, singleQuality);
\r
240 // if QueryBases is empty (and this is a allowed case)
\r
241 if ( bAlignment.QueryBases.empty() )
\r
242 bAlignment.AlignedBases = bAlignment.QueryBases;
\r
244 // if QueryBases contains data, then build AlignedBases using CIGAR data
\r
247 // resize AlignedBases
\r
248 bAlignment.AlignedBases.clear();
\r
249 bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);
\r
251 // iterate over CigarOps
\r
253 vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();
\r
254 vector<CigarOp>::const_iterator cigarEnd = bAlignment.CigarData.end();
\r
255 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
\r
257 const CigarOp& op = (*cigarIter);
\r
262 bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases
\r
266 k += op.Length; // for 'S' - soft clip, skip over query bases
\r
270 bAlignment.AlignedBases.append(op.Length, '-'); // for 'D' - write gap character
\r
274 bAlignment.AlignedBases.append( op.Length, '*' ); // for 'P' - write padding character
\r
278 bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in original query sequence
\r
282 break; // for 'H' - hard clip, do nothing to AlignedBases, move to next op
\r
285 printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here
\r
291 // -----------------------
\r
292 // Added: 3-25-2010 DB
\r
293 // Fixed: endian-correctness for tag data
\r
294 // -----------------------
\r
295 if ( IsBigEndian ) {
\r
297 while ( (unsigned int)i < tagDataLength ) {
\r
299 i += 2; // skip tag type (e.g. "RG", "NM", etc)
\r
300 uint8_t type = toupper(tagData[i]); // lower & upper case letters have same meaning
\r
301 ++i; // skip value type
\r
311 SwapEndian_16p(&tagData[i]);
\r
312 i += sizeof(uint16_t);
\r
317 SwapEndian_32p(&tagData[i]);
\r
318 i += sizeof(uint32_t);
\r
322 SwapEndian_64p(&tagData[i]);
\r
323 i += sizeof(uint64_t);
\r
328 while (tagData[i]) { ++i; }
\r
329 ++i; // increment one more for null terminator
\r
333 printf("ERROR: Invalid tag value type\n"); // shouldn't get here
\r
340 bAlignment.TagData.clear();
\r
341 bAlignment.TagData.resize(tagDataLength);
\r
342 memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);
\r
344 // clear the core-only flag
\r
345 bAlignment.SupportData.HasCoreOnly = false;
\r
351 // clear index data structure
\r
352 void BamReader::BamReaderPrivate::ClearIndex(void) {
\r
355 IsIndexLoaded = false;
\r
358 // closes the BAM file
\r
359 void BamReader::BamReaderPrivate::Close(void) {
\r
361 // close BGZF file stream
\r
364 // clear out index data
\r
367 // clear out header data
\r
368 HeaderText.clear();
\r
370 // clear out region flags
\r
371 IsLeftBoundSpecified = false;
\r
372 IsRightBoundSpecified = false;
\r
373 IsRegionSpecified = false;
\r
376 // creates index for BAM file, saves to file
\r
377 // default behavior is to create the BAM standard index (".bai")
\r
378 // set flag to false to create the BamTools-specific index (".bti")
\r
379 bool BamReader::BamReaderPrivate::CreateIndex(bool useStandardIndex) {
\r
381 // clear out prior index data
\r
384 // create index based on type requested
\r
385 if ( useStandardIndex )
\r
386 NewIndex = new BamDefaultIndex(&mBGZF, Parent, IsBigEndian);
\r
387 // create BamTools 'custom' index
\r
389 NewIndex = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);
\r
393 ok &= NewIndex->Build();
\r
394 IsIndexLoaded = ok;
\r
396 // attempt to save index data to file
\r
397 ok &= NewIndex->Write(Filename);
\r
399 // return success/fail of both building & writing index
\r
403 // get next alignment (from specified region, if given)
\r
404 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {
\r
406 // if valid alignment found, attempt to parse char data, and return success/failure
\r
407 if ( GetNextAlignmentCore(bAlignment) )
\r
408 return BuildCharData(bAlignment);
\r
410 // no valid alignment found
\r
415 // retrieves next available alignment core data (returns success/fail)
\r
416 // ** DOES NOT parse any character data (read name, bases, qualities, tag data)
\r
417 // these can be accessed, if necessary, from the supportData
\r
418 // useful for operations requiring ONLY positional or other alignment-related information
\r
419 bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {
\r
421 // if valid alignment available
\r
422 if ( LoadNextAlignment(bAlignment) ) {
\r
424 // set core-only flag
\r
425 bAlignment.SupportData.HasCoreOnly = true;
\r
427 // if region not specified, return success
\r
428 if ( !IsLeftBoundSpecified ) return true;
\r
430 // determine region state (before, within, after)
\r
431 BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);
\r
433 // if alignment lies after region, return false
\r
434 if ( state == AFTER_REGION ) return false;
\r
436 while ( state != WITHIN_REGION ) {
\r
437 // if no valid alignment available (likely EOF) return failure
\r
438 if ( !LoadNextAlignment(bAlignment) ) return false;
\r
439 // if alignment lies after region, return false (no available read within region)
\r
440 state = IsOverlap(bAlignment);
\r
441 if ( state == AFTER_REGION ) return false;
\r
444 // return success (alignment found that overlaps region)
\r
448 // no valid alignment
\r
453 // returns RefID for given RefName (returns References.size() if not found)
\r
454 int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {
\r
456 // retrieve names from reference data
\r
457 vector<string> refNames;
\r
458 RefVector::const_iterator refIter = References.begin();
\r
459 RefVector::const_iterator refEnd = References.end();
\r
460 for ( ; refIter != refEnd; ++refIter)
\r
461 refNames.push_back( (*refIter).RefName );
\r
463 // return 'index-of' refName ( if not found, returns refNames.size() )
\r
464 return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
\r
467 // returns region state - whether alignment ends before, overlaps, or starts after currently specified region
\r
468 // this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
\r
469 BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {
\r
471 // --------------------------------------------------
\r
472 // check alignment start against right bound cutoff
\r
474 // if full region of interest was given
\r
475 if ( IsRightBoundSpecified ) {
\r
477 // read starts on right bound reference, but AFTER right bound position
\r
478 if ( bAlignment.RefID == Region.RightRefID && bAlignment.Position > Region.RightPosition )
\r
479 return AFTER_REGION;
\r
481 // if read starts on reference AFTER right bound, return false
\r
482 if ( bAlignment.RefID > Region.RightRefID )
\r
483 return AFTER_REGION;
\r
486 // --------------------------------------------------------
\r
487 // no right bound given OR read starts before right bound
\r
488 // so, check if it overlaps left bound
\r
490 // if read starts on left bound reference AND after left boundary, return success
\r
491 if ( bAlignment.RefID == Region.LeftRefID && bAlignment.Position >= Region.LeftPosition)
\r
492 return WITHIN_REGION;
\r
494 // if read is on any reference sequence before left bound, return false
\r
495 if ( bAlignment.RefID < Region.LeftRefID )
\r
496 return BEFORE_REGION;
\r
498 // --------------------------------------------------------
\r
499 // read is on left bound reference, but starts before left bound position
\r
501 // if it overlaps, return WITHIN_REGION
\r
502 if ( bAlignment.GetEndPosition() >= Region.LeftPosition )
\r
503 return WITHIN_REGION;
\r
504 // else begins before left bound position
\r
506 return BEFORE_REGION;
\r
509 // jumps to specified region(refID, leftBound) in BAM file, returns success/fail
\r
510 bool BamReader::BamReaderPrivate::Jump(int refID, int position) {
\r
512 // -----------------------------------------------------------------------
\r
513 // check for existing index
\r
514 if ( !IsIndexLoaded || NewIndex == 0 ) return false;
\r
515 // see if reference has alignments
\r
516 if ( !NewIndex->HasAlignments(refID) ) return false;
\r
517 // make sure position is valid
\r
518 if ( position > References.at(refID).RefLength ) return false;
\r
520 // determine possible offsets
\r
521 vector<int64_t> offsets;
\r
522 if ( !NewIndex->GetOffsets(Region, IsRightBoundSpecified, offsets) ) {
\r
523 printf("ERROR: Could not jump: unable to calculate offset for specified region.\n");
\r
527 // iterate through offsets
\r
528 BamAlignment bAlignment;
\r
529 bool result = true;
\r
530 for ( vector<int64_t>::const_iterator o = offsets.begin(); o != offsets.end(); ++o) {
\r
532 // attempt seek & load first available alignment
\r
533 result &= mBGZF.Seek(*o);
\r
534 LoadNextAlignment(bAlignment);
\r
536 // if this alignment corresponds to desired position
\r
537 // return success of seeking back to 'current offset'
\r
538 if ( (bAlignment.RefID == refID && bAlignment.Position + bAlignment.Length > position) || (bAlignment.RefID > refID) ) {
\r
539 if ( o != offsets.begin() ) --o;
\r
540 return mBGZF.Seek(*o);
\r
547 // load BAM header data
\r
548 void BamReader::BamReaderPrivate::LoadHeaderData(void) {
\r
550 // check to see if proper BAM header
\r
552 if (mBGZF.Read(buffer, 4) != 4) {
\r
553 printf("Could not read header type\n");
\r
557 if (strncmp(buffer, "BAM\001", 4)) {
\r
558 printf("wrong header type!\n");
\r
562 // get BAM header text length
\r
563 mBGZF.Read(buffer, 4);
\r
564 unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);
\r
565 if ( IsBigEndian ) SwapEndian_32(headerTextLength);
\r
567 // get BAM header text
\r
568 char* headerText = (char*)calloc(headerTextLength + 1, 1);
\r
569 mBGZF.Read(headerText, headerTextLength);
\r
570 HeaderText = (string)((const char*)headerText);
\r
572 // clean up calloc-ed temp variable
\r
576 // load existing index data from BAM index file (".bti" OR ".bai"), return success/fail
\r
577 bool BamReader::BamReaderPrivate::LoadIndex(const bool lookForIndex, const bool preferStandardIndex) {
\r
579 // clear out any existing index data
\r
582 // if no index filename provided, so we need to look for available index files
\r
583 if ( IndexFilename.empty() ) {
\r
585 // attempt to load BamIndex based on current Filename provided & preferStandardIndex flag
\r
586 const BamIndex::PreferredIndexType type = (preferStandardIndex ? BamIndex::STANDARD : BamIndex::BAMTOOLS);
\r
587 NewIndex = BamIndex::FromBamFilename(Filename, &mBGZF, Parent, IsBigEndian, type);
\r
589 // if null, return failure
\r
590 if ( NewIndex == 0 ) return false;
\r
592 // generate proper IndexFilename based on type of index created
\r
593 IndexFilename = Filename + NewIndex->Extension();
\r
597 // attempt to load BamIndex based on IndexFilename provided by client
\r
598 NewIndex = BamIndex::FromIndexFilename(IndexFilename, &mBGZF, Parent, IsBigEndian);
\r
600 // if null, return failure
\r
601 if ( NewIndex == 0 ) return false;
\r
604 // an index file was found
\r
605 // return success of loading the index data from file
\r
606 IsIndexLoaded = NewIndex->Load(IndexFilename);
\r
607 return IsIndexLoaded;
\r
610 // populates BamAlignment with alignment data under file pointer, returns success/fail
\r
611 bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {
\r
613 // read in the 'block length' value, make sure it's not zero
\r
615 mBGZF.Read(buffer, 4);
\r
616 bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);
\r
617 if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }
\r
618 if ( bAlignment.SupportData.BlockLength == 0 ) return false;
\r
620 // read in core alignment data, make sure the right size of data was read
\r
621 char x[BAM_CORE_SIZE];
\r
622 if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) return false;
\r
624 if ( IsBigEndian ) {
\r
625 for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) )
\r
626 SwapEndian_32p(&x[i]);
\r
629 // set BamAlignment 'core' and 'support' data
\r
630 bAlignment.RefID = BgzfData::UnpackSignedInt(&x[0]);
\r
631 bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);
\r
633 unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);
\r
634 bAlignment.Bin = tempValue >> 16;
\r
635 bAlignment.MapQuality = tempValue >> 8 & 0xff;
\r
636 bAlignment.SupportData.QueryNameLength = tempValue & 0xff;
\r
638 tempValue = BgzfData::UnpackUnsignedInt(&x[12]);
\r
639 bAlignment.AlignmentFlag = tempValue >> 16;
\r
640 bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;
\r
642 bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);
\r
643 bAlignment.MateRefID = BgzfData::UnpackSignedInt(&x[20]);
\r
644 bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);
\r
645 bAlignment.InsertSize = BgzfData::UnpackSignedInt(&x[28]);
\r
647 // set BamAlignment length
\r
648 bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;
\r
650 // read in character data - make sure proper data size was read
\r
651 bool readCharDataOK = false;
\r
652 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
\r
653 char* allCharData = (char*)calloc(sizeof(char), dataLength);
\r
655 if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) {
\r
657 // store 'allCharData' in supportData structure
\r
658 bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
\r
660 // set success flag
\r
661 readCharDataOK = true;
\r
664 // need to calculate this here so that BamAlignment::GetEndPosition() performs correctly,
\r
665 // even when BamReader::GetNextAlignmentCore() is called
\r
666 const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;
\r
667 uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
\r
669 bAlignment.CigarData.clear();
\r
670 bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);
\r
671 for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {
\r
673 // swap if necessary
\r
674 if ( IsBigEndian ) SwapEndian_32(cigarData[i]);
\r
676 // build CigarOp structure
\r
677 op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
\r
678 op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
\r
681 bAlignment.CigarData.push_back(op);
\r
686 return readCharDataOK;
\r
689 // loads reference data from BAM file
\r
690 void BamReader::BamReaderPrivate::LoadReferenceData(void) {
\r
692 // get number of reference sequences
\r
694 mBGZF.Read(buffer, 4);
\r
695 unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);
\r
696 if ( IsBigEndian ) SwapEndian_32(numberRefSeqs);
\r
697 if ( numberRefSeqs == 0 ) return;
\r
698 References.reserve((int)numberRefSeqs);
\r
700 // iterate over all references in header
\r
701 for (unsigned int i = 0; i != numberRefSeqs; ++i) {
\r
703 // get length of reference name
\r
704 mBGZF.Read(buffer, 4);
\r
705 unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);
\r
706 if ( IsBigEndian ) SwapEndian_32(refNameLength);
\r
707 char* refName = (char*)calloc(refNameLength, 1);
\r
709 // get reference name and reference sequence length
\r
710 mBGZF.Read(refName, refNameLength);
\r
711 mBGZF.Read(buffer, 4);
\r
712 int refLength = BgzfData::UnpackSignedInt(buffer);
\r
713 if ( IsBigEndian ) SwapEndian_32(refLength);
\r
715 // store data for reference
\r
716 RefData aReference;
\r
717 aReference.RefName = (string)((const char*)refName);
\r
718 aReference.RefLength = refLength;
\r
719 References.push_back(aReference);
\r
721 // clean up calloc-ed temp variable
\r
726 // opens BAM file (and index)
\r
727 bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename, const bool lookForIndex, const bool preferStandardIndex) {
\r
730 Filename = filename;
\r
731 IndexFilename = indexFilename;
\r
733 // open the BGZF file for reading, return false on failure
\r
734 if ( !mBGZF.Open(filename, "rb") ) return false;
\r
736 // retrieve header text & reference data
\r
738 LoadReferenceData();
\r
740 // store file offset of first alignment
\r
741 AlignmentsBeginOffset = mBGZF.Tell();
\r
743 // if no index filename provided
\r
744 if ( IndexFilename.empty() ) {
\r
746 // client did not specify that index SHOULD be found
\r
747 // useful for cases where sequential access is all that is required
\r
748 if ( !lookForIndex ) return true;
\r
750 // otherwise, look for index file, return success/fail
\r
751 return LoadIndex(lookForIndex, preferStandardIndex) ;
\r
754 // client supplied an index filename
\r
755 // attempt to load index data, return success/fail
\r
756 return LoadIndex(lookForIndex, preferStandardIndex);
\r
759 // returns BAM file pointer to beginning of alignment data
\r
760 bool BamReader::BamReaderPrivate::Rewind(void) {
\r
762 // rewind to first alignment
\r
763 if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;
\r
765 // retrieve first alignment data
\r
767 if ( !LoadNextAlignment(al) ) return false;
\r
769 // reset default region info using first alignment in file
\r
770 Region.LeftRefID = al.RefID;
\r
771 Region.LeftPosition = al.Position;
\r
772 Region.RightRefID = -1;
\r
773 Region.RightPosition = -1;
\r
774 IsLeftBoundSpecified = false;
\r
775 IsRightBoundSpecified = false;
\r
777 // rewind back to before first alignment
\r
778 // return success/fail of seek
\r
779 return mBGZF.Seek(AlignmentsBeginOffset);
\r
782 // sets a region of interest (with left & right bound reference/position)
\r
783 // attempts a Jump() to left bound as well
\r
784 // returns success/failure of Jump()
\r
785 bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {
\r
787 // save region of interest
\r
791 if ( region.LeftRefID >= 0 && region.LeftPosition >= 0 ) IsLeftBoundSpecified = true;
\r
792 if ( region.RightRefID >= 0 && region.RightPosition >= 0 ) IsRightBoundSpecified = true;
\r
794 // attempt jump to beginning of region, return success/fail of Jump()
\r
795 return Jump( Region.LeftRefID, Region.LeftPosition );
\r