1 // ***************************************************************************
2 // BamReader_p.cpp (c) 2009 Derek Barnett
3 // Marth Lab, Department of Biology, Boston College
4 // All rights reserved.
5 // ---------------------------------------------------------------------------
6 // Last modified: 22 November 2010 (DB)
7 // ---------------------------------------------------------------------------
8 // Provides the basic functionality for reading BAM files
9 // ***************************************************************************
11 #include <api/BamReader.h>
13 #include <api/internal/BamReader_p.h>
14 #include <api/internal/BamStandardIndex_p.h>
15 #include <api/internal/BamToolsIndex_p.h>
16 using namespace BamTools;
17 using namespace BamTools::Internal;
26 BamReaderPrivate::BamReaderPrivate(BamReader* parent)
30 , AlignmentsBeginOffset(0)
32 , IndexCacheMode(BamIndex::LimitedIndexCaching)
33 , HasAlignmentsInRegion(true)
35 , DNA_LOOKUP("=ACMGRSVTWYHKDBN")
36 , CIGAR_LOOKUP("MIDNSHP")
38 IsBigEndian = SystemIsBigEndian();
42 BamReaderPrivate::~BamReaderPrivate(void) {
46 // adjusts requested region if necessary (depending on where data actually begins)
47 void BamReaderPrivate::AdjustRegion(BamRegion& region) {
49 // check for valid index first
50 if ( Index == 0 ) return;
52 // see if any references in region have alignments
53 HasAlignmentsInRegion = false;
54 int currentId = region.LeftRefID;
56 const int rightBoundRefId = ( region.isRightBoundSpecified() ? region.RightRefID : References.size() - 1 );
57 while ( currentId <= rightBoundRefId ) {
58 HasAlignmentsInRegion = Index->HasAlignments(currentId);
59 if ( HasAlignmentsInRegion ) break;
63 // if no data found on any reference in region
64 if ( !HasAlignmentsInRegion ) return;
66 // if left bound of desired region had no data, use first reference that had data
67 // otherwise, leave requested region as-is
68 if ( currentId != region.LeftRefID ) {
69 region.LeftRefID = currentId;
70 region.LeftPosition = 0;
74 // fills out character data for BamAlignment data
75 bool BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {
77 // calculate character lengths/offsets
78 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
79 const unsigned int seqDataOffset = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);
80 const unsigned int qualDataOffset = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;
81 const unsigned int tagDataOffset = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;
82 const unsigned int tagDataLength = dataLength - tagDataOffset;
84 // set up char buffers
85 const char* allCharData = bAlignment.SupportData.AllCharData.data();
86 const char* seqData = ((const char*)allCharData) + seqDataOffset;
87 const char* qualData = ((const char*)allCharData) + qualDataOffset;
88 char* tagData = ((char*)allCharData) + tagDataOffset;
90 // store alignment name (relies on null char in name as terminator)
91 bAlignment.Name.assign((const char*)(allCharData));
93 // save query sequence
94 bAlignment.QueryBases.clear();
95 bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);
96 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
97 char singleBase = DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
98 bAlignment.QueryBases.append(1, singleBase);
101 // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
102 bAlignment.Qualities.clear();
103 bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);
104 for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
105 char singleQuality = (char)(qualData[i]+33);
106 bAlignment.Qualities.append(1, singleQuality);
109 // if QueryBases is empty (and this is a allowed case)
110 if ( bAlignment.QueryBases.empty() )
111 bAlignment.AlignedBases = bAlignment.QueryBases;
113 // if QueryBases contains data, then build AlignedBases using CIGAR data
116 // resize AlignedBases
117 bAlignment.AlignedBases.clear();
118 bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);
120 // iterate over CigarOps
122 vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();
123 vector<CigarOp>::const_iterator cigarEnd = bAlignment.CigarData.end();
124 for ( ; cigarIter != cigarEnd; ++cigarIter ) {
126 const CigarOp& op = (*cigarIter);
131 bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases
135 k += op.Length; // for 'S' - soft clip, skip over query bases
139 bAlignment.AlignedBases.append(op.Length, '-'); // for 'D' - write gap character
143 bAlignment.AlignedBases.append( op.Length, '*' ); // for 'P' - write padding character
147 bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in original query sequence
151 break; // for 'H' - hard clip, do nothing to AlignedBases, move to next op
154 fprintf(stderr, "ERROR: Invalid Cigar op type\n"); // shouldn't get here
160 // -----------------------
161 // Added: 3-25-2010 DB
162 // Fixed: endian-correctness for tag data
163 // -----------------------
166 while ( (unsigned int)i < tagDataLength ) {
168 i += 2; // skip tag type (e.g. "RG", "NM", etc)
169 uint8_t type = toupper(tagData[i]); // lower & upper case letters have same meaning
170 ++i; // skip value type
180 SwapEndian_16p(&tagData[i]);
181 i += sizeof(uint16_t);
186 SwapEndian_32p(&tagData[i]);
187 i += sizeof(uint32_t);
191 SwapEndian_64p(&tagData[i]);
192 i += sizeof(uint64_t);
197 while (tagData[i]) { ++i; }
198 ++i; // increment one more for null terminator
202 fprintf(stderr, "ERROR: Invalid tag value type\n"); // shouldn't get here
209 bAlignment.TagData.clear();
210 bAlignment.TagData.resize(tagDataLength);
211 memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);
213 // clear the core-only flag
214 bAlignment.SupportData.HasCoreOnly = false;
220 // clear index data structure
221 void BamReaderPrivate::ClearIndex(void) {
227 // closes the BAM file
228 void BamReaderPrivate::Close(void) {
230 // close BGZF file stream
233 // clear out index data
236 // clear out header data
243 // clear out region flags
247 // creates index for BAM file, saves to file
248 // default behavior is to create the BAM standard index (".bai")
249 // set flag to false to create the BamTools-specific index (".bti")
250 bool BamReaderPrivate::CreateIndex(bool useStandardIndex) {
252 // clear out prior index data
255 // create index based on type requested
256 if ( useStandardIndex )
257 Index = new BamStandardIndex(&mBGZF, Parent);
259 Index = new BamToolsIndex(&mBGZF, Parent);
261 // set index cache mode to full for writing
262 Index->SetCacheMode(BamIndex::FullIndexCaching);
266 ok &= Index->Build();
269 // mark empty references
272 // attempt to save index data to file
273 ok &= Index->Write(Filename);
275 // set client's desired index cache mode
276 Index->SetCacheMode(IndexCacheMode);
278 // return success/fail of both building & writing index
282 const string BamReaderPrivate::GetHeaderText(void) const {
287 // return m_header->Text();
289 // return string("");
292 // get next alignment (from specified region, if given)
293 bool BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {
295 // if valid alignment found, attempt to parse char data, and return success/failure
296 if ( GetNextAlignmentCore(bAlignment) )
297 return BuildCharData(bAlignment);
299 // no valid alignment found
303 // retrieves next available alignment core data (returns success/fail)
304 // ** DOES NOT parse any character data (read name, bases, qualities, tag data)
305 // these can be accessed, if necessary, from the supportData
306 // useful for operations requiring ONLY positional or other alignment-related information
307 bool BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {
309 // if region is set but has no alignments
310 if ( !Region.isNull() && !HasAlignmentsInRegion )
313 // if valid alignment available
314 if ( LoadNextAlignment(bAlignment) ) {
316 // set core-only flag
317 bAlignment.SupportData.HasCoreOnly = true;
319 // if region not specified with at least a left boundary, return success
320 if ( !Region.isLeftBoundSpecified() ) return true;
322 // determine region state (before, within, after)
323 BamReaderPrivate::RegionState state = IsOverlap(bAlignment);
325 // if alignment lies after region, return false
326 if ( state == AFTER_REGION ) return false;
328 while ( state != WITHIN_REGION ) {
329 // if no valid alignment available (likely EOF) return failure
330 if ( !LoadNextAlignment(bAlignment) ) return false;
331 // if alignment lies after region, return false (no available read within region)
332 state = IsOverlap(bAlignment);
333 if ( state == AFTER_REGION ) return false;
336 // return success (alignment found that overlaps region)
340 // no valid alignment
344 // returns RefID for given RefName (returns References.size() if not found)
345 int BamReaderPrivate::GetReferenceID(const string& refName) const {
347 // retrieve names from reference data
348 vector<string> refNames;
349 RefVector::const_iterator refIter = References.begin();
350 RefVector::const_iterator refEnd = References.end();
351 for ( ; refIter != refEnd; ++refIter)
352 refNames.push_back( (*refIter).RefName );
354 // return 'index-of' refName ( if not found, returns refNames.size() )
355 return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
358 // returns region state - whether alignment ends before, overlaps, or starts after currently specified region
359 // this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
360 BamReaderPrivate::RegionState BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {
362 // if alignment is on any reference sequence before left bound
363 if ( bAlignment.RefID < Region.LeftRefID ) return BEFORE_REGION;
365 // if alignment starts on left bound reference
366 else if ( bAlignment.RefID == Region.LeftRefID ) {
368 // if alignment starts at or after left boundary
369 if ( bAlignment.Position >= Region.LeftPosition) {
371 // if right boundary is specified AND
372 // left/right boundaries are on same reference AND
373 // alignment starts past right boundary
374 if ( Region.isRightBoundSpecified() &&
375 Region.LeftRefID == Region.RightRefID &&
376 bAlignment.Position > Region.RightPosition )
379 // otherwise, alignment is within region
380 return WITHIN_REGION;
383 // alignment starts before left boundary
385 // check if alignment overlaps left boundary
386 if ( bAlignment.GetEndPosition() >= Region.LeftPosition ) return WITHIN_REGION;
387 else return BEFORE_REGION;
391 // alignment starts on a reference after the left bound
394 // if region has a right boundary
395 if ( Region.isRightBoundSpecified() ) {
397 // alignment is on reference between boundaries
398 if ( bAlignment.RefID < Region.RightRefID ) return WITHIN_REGION;
400 // alignment is on reference after right boundary
401 else if ( bAlignment.RefID > Region.RightRefID ) return AFTER_REGION;
403 // alignment is on right bound reference
405 // check if alignment starts before or at right boundary
406 if ( bAlignment.Position <= Region.RightPosition ) return WITHIN_REGION;
407 else return AFTER_REGION;
411 // otherwise, alignment is after left bound reference, but there is no right boundary
412 else return WITHIN_REGION;
416 // load BAM header data
417 void BamReaderPrivate::LoadHeaderData(void) {
419 // m_header = new BamHeader(&mBGZF);
420 // bool headerLoadedOk = m_header->Load();
421 // if ( !headerLoadedOk )
422 // cerr << "BamReader could not load header" << endl;
424 // check to see if proper BAM header
426 if (mBGZF.Read(buffer, 4) != 4) {
427 fprintf(stderr, "Could not read header type\n");
431 if (strncmp(buffer, "BAM\001", 4)) {
432 fprintf(stderr, "wrong header type!\n");
436 // get BAM header text length
437 mBGZF.Read(buffer, 4);
438 unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);
439 if ( IsBigEndian ) SwapEndian_32(headerTextLength);
441 // get BAM header text
442 char* headerText = (char*)calloc(headerTextLength + 1, 1);
443 mBGZF.Read(headerText, headerTextLength);
444 HeaderText = (string)((const char*)headerText);
446 // clean up calloc-ed temp variable
450 // load existing index data from BAM index file (".bti" OR ".bai"), return success/fail
451 bool BamReaderPrivate::LoadIndex(const bool lookForIndex, const bool preferStandardIndex) {
453 // clear out any existing index data
456 // if no index filename provided, so we need to look for available index files
457 if ( IndexFilename.empty() ) {
459 // attempt to load BamIndex based on current Filename provided & preferStandardIndex flag
460 const BamIndex::PreferredIndexType type = (preferStandardIndex ? BamIndex::STANDARD : BamIndex::BAMTOOLS);
461 Index = BamIndex::FromBamFilename(Filename, &mBGZF, Parent, type);
463 // if null, return failure
464 if ( Index == 0 ) return false;
466 // generate proper IndexFilename based on type of index created
467 IndexFilename = Filename + Index->Extension();
472 // attempt to load BamIndex based on IndexFilename provided by client
473 Index = BamIndex::FromIndexFilename(IndexFilename, &mBGZF, Parent);
475 // if null, return failure
476 if ( Index == 0 ) return false;
479 // set cache mode for BamIndex
480 Index->SetCacheMode(IndexCacheMode);
482 // loading the index data from file
483 HasIndex = Index->Load(IndexFilename);
485 // mark empty references
488 // return index status
492 // populates BamAlignment with alignment data under file pointer, returns success/fail
493 bool BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {
495 // read in the 'block length' value, make sure it's not zero
497 mBGZF.Read(buffer, 4);
498 bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);
499 if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }
500 if ( bAlignment.SupportData.BlockLength == 0 ) return false;
502 // read in core alignment data, make sure the right size of data was read
503 char x[BAM_CORE_SIZE];
504 if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) return false;
507 for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) )
508 SwapEndian_32p(&x[i]);
511 // set BamAlignment 'core' and 'support' data
512 bAlignment.RefID = BgzfData::UnpackSignedInt(&x[0]);
513 bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);
515 unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);
516 bAlignment.Bin = tempValue >> 16;
517 bAlignment.MapQuality = tempValue >> 8 & 0xff;
518 bAlignment.SupportData.QueryNameLength = tempValue & 0xff;
520 tempValue = BgzfData::UnpackUnsignedInt(&x[12]);
521 bAlignment.AlignmentFlag = tempValue >> 16;
522 bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;
524 bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);
525 bAlignment.MateRefID = BgzfData::UnpackSignedInt(&x[20]);
526 bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);
527 bAlignment.InsertSize = BgzfData::UnpackSignedInt(&x[28]);
529 // set BamAlignment length
530 bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;
532 // read in character data - make sure proper data size was read
533 bool readCharDataOK = false;
534 const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
535 char* allCharData = (char*)calloc(sizeof(char), dataLength);
537 if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) {
539 // store 'allCharData' in supportData structure
540 bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
543 readCharDataOK = true;
546 // need to calculate this here so that BamAlignment::GetEndPosition() performs correctly,
547 // even when GetNextAlignmentCore() is called
548 const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;
549 uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
551 bAlignment.CigarData.clear();
552 bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);
553 for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {
556 if ( IsBigEndian ) SwapEndian_32(cigarData[i]);
558 // build CigarOp structure
559 op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
560 op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
563 bAlignment.CigarData.push_back(op);
568 return readCharDataOK;
571 // loads reference data from BAM file
572 void BamReaderPrivate::LoadReferenceData(void) {
574 // get number of reference sequences
576 mBGZF.Read(buffer, 4);
577 unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);
578 if ( IsBigEndian ) SwapEndian_32(numberRefSeqs);
579 if ( numberRefSeqs == 0 ) return;
580 References.reserve((int)numberRefSeqs);
582 // iterate over all references in header
583 for (unsigned int i = 0; i != numberRefSeqs; ++i) {
585 // get length of reference name
586 mBGZF.Read(buffer, 4);
587 unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);
588 if ( IsBigEndian ) SwapEndian_32(refNameLength);
589 char* refName = (char*)calloc(refNameLength, 1);
591 // get reference name and reference sequence length
592 mBGZF.Read(refName, refNameLength);
593 mBGZF.Read(buffer, 4);
594 int refLength = BgzfData::UnpackSignedInt(buffer);
595 if ( IsBigEndian ) SwapEndian_32(refLength);
597 // store data for reference
599 aReference.RefName = (string)((const char*)refName);
600 aReference.RefLength = refLength;
601 References.push_back(aReference);
603 // clean up calloc-ed temp variable
608 // mark references with no alignment data
609 void BamReaderPrivate::MarkReferences(void) {
611 // ensure index is available
612 if ( !HasIndex ) return;
614 // mark empty references
615 for ( int i = 0; i < (int)References.size(); ++i )
616 References.at(i).RefHasAlignments = Index->HasAlignments(i);
619 // opens BAM file (and index)
620 bool BamReaderPrivate::Open(const string& filename, const string& indexFilename, const bool lookForIndex, const bool preferStandardIndex) {
624 IndexFilename = indexFilename;
626 // open the BGZF file for reading, return false on failure
627 if ( !mBGZF.Open(filename, "rb") ) return false;
629 // retrieve header text & reference data
633 // store file offset of first alignment
634 AlignmentsBeginOffset = mBGZF.Tell();
636 // if no index filename provided
637 if ( IndexFilename.empty() ) {
639 // client did not specify that index SHOULD be found
640 // useful for cases where sequential access is all that is required
641 if ( !lookForIndex ) return true;
643 // otherwise, look for index file, return success/fail
644 return LoadIndex(lookForIndex, preferStandardIndex) ;
647 // client supplied an index filename
648 // attempt to load index data, return success/fail
649 return LoadIndex(lookForIndex, preferStandardIndex);
652 // returns BAM file pointer to beginning of alignment data
653 bool BamReaderPrivate::Rewind(void) {
655 // rewind to first alignment, return false if unable to seek
656 if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;
658 // retrieve first alignment data, return false if unable to read
660 if ( !LoadNextAlignment(al) ) return false;
662 // reset default region info using first alignment in file
664 HasAlignmentsInRegion = true;
666 // rewind back to beginning of first alignment
667 // return success/fail of seek
668 return mBGZF.Seek(AlignmentsBeginOffset);
671 // change the index caching behavior
672 void BamReaderPrivate::SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode) {
673 IndexCacheMode = mode;
674 if ( Index == 0 ) return;
675 Index->SetCacheMode(mode);
678 // asks Index to attempt a Jump() to specified region
679 // returns success/failure
680 bool BamReaderPrivate::SetRegion(const BamRegion& region) {
682 // clear out any prior BamReader region data
684 // N.B. - this is cleared so that BamIndex now has free reign to call
685 // GetNextAlignmentCore() and do overlap checking without worrying about BamReader
686 // performing any overlap checking of its own and moving on to the next read... Calls
687 // to GetNextAlignmentCore() with no Region set, simply return the next alignment.
688 // This ensures that the Index is able to do just that. (All without exposing
689 // LoadNextAlignment() to the public API, and potentially confusing clients with the nomenclature)
692 // check for existing index
693 if ( !HasIndex ) return false;
695 // adjust region if necessary to reflect where data actually begins
696 BamRegion adjustedRegion(region);
697 AdjustRegion(adjustedRegion);
699 // if no data present, return true
700 // not an error, but BamReader knows that no data is there for future alignment access
701 // (this is useful in a MultiBamReader setting where some BAM files may lack data in regions
702 // that other BAMs have data)
703 if ( !HasAlignmentsInRegion ) {
704 Region = adjustedRegion;
708 // attempt jump to user-specified region return false if jump could not be performed at all
709 // (invalid index, unknown reference, etc)
711 // Index::Jump() is allowed to modify the HasAlignmentsInRegion flag
712 // * This covers case where a region is requested that lies beyond the last alignment on a reference
713 // If this occurs, any subsequent calls to GetNexAlignment[Core] simply return false
714 // BamMultiReader is then able to successfully pull alignments from a region from multiple files
715 // even if one or more have no data.
716 if ( !Index->Jump(adjustedRegion, &HasAlignmentsInRegion) ) return false;
718 // save region and return success
719 Region = adjustedRegion;