4 // Marth Lab, Boston College
5 // Last modified: 23 April 2009
12 // static character constants
13 const char* BamReader::DNA_LOOKUP = "=ACMGRSVTWYHKDBN";
14 const char* BamReader::CIGAR_LOOKUP = "MIDNSHP";
16 BamReader::BamReader(const char* filename, const char* indexFilename)
17 : m_filename( (char*)filename )
18 , m_indexFilename( (char*)indexFilename )
23 , m_isIndexLoaded(false)
24 , m_isRegionSpecified(false)
25 , m_isCalculateAlignedBases(true)
28 , m_alignmentsBeginOffset(0)
33 BamReader::~BamReader(void) {
36 if ( m_isOpen ) { Close(); }
40 bool BamReader::Open(void) {
42 if (!m_isOpen && m_filename != NULL ) {
45 m_file = bam_open(m_filename, "r");
47 // get header info && index info
48 if ( (m_file != NULL) && LoadHeader() ) {
50 // save file offset where alignments start
51 m_alignmentsBeginOffset = bam_tell(m_file);
57 // try to open (and load) index data, if index file given
58 if ( m_indexFilename != NULL ) {
66 bool BamReader::OpenIndex(void) {
68 if ( m_indexFilename && !m_isIndexLoaded ) {
69 m_isIndexLoaded = LoadIndex();
71 return m_isIndexLoaded;
75 bool BamReader::Close(void) {
80 int ret = bam_close(m_file);
83 if ( m_index != NULL) { delete m_index; }
89 m_isIndexLoaded = false;
92 m_isRegionSpecified = false;
94 // return success/fail of bam_close
102 const char* BamReader::Filename(void) const {
103 return (const char*)m_filename;
107 void BamReader::SetFilename(const char* filename) {
108 m_filename = (char*)filename;
111 // get BAM Index filename
112 const char* BamReader::IndexFilename(void) const {
113 return (const char*)m_indexFilename;
116 // set BAM Index filename
117 void BamReader::SetIndexFilename(const char* indexFilename) {
118 m_indexFilename = (char*)indexFilename;
121 // return full header text
122 const string BamReader::GetHeaderText(void) const {
126 // return number of reference sequences in BAM file
127 const int BamReader::GetReferenceCount(void) const {
128 return m_references.size();
131 // get RefID from reference name
132 const int BamReader::GetRefID(string refName) const {
134 vector<string> refNames;
135 RefVector::const_iterator refIter = m_references.begin();
136 RefVector::const_iterator refEnd = m_references.end();
137 for ( ; refIter != refEnd; ++refIter) {
138 refNames.push_back( (*refIter).RefName );
141 // return 'index-of' refName (if not found, returns refNames.size())
142 return Index( refNames.begin(), refNames.end(), refName );
145 const RefVector BamReader::GetReferenceData(void) const {
149 bool BamReader::Jump(int refID, unsigned int left) {
151 // if index available, and region is valid
152 if ( m_isIndexLoaded && m_references.at(refID).RefHasAlignments && (left <= m_references.at(refID).RefLength) ) {
153 m_currentRefID = refID;
154 m_currentLeft = left;
155 m_isRegionSpecified = true;
157 int64_t offset = GetOffset(m_currentRefID, m_currentLeft);
158 if ( offset == -1 ) { return false; }
159 else { return ( bam_seek(m_file, offset, SEEK_SET) == 0 ); }
164 bool BamReader::Rewind(void) {
167 int refCount = m_references.size();
168 for ( ; refID < refCount; ++refID ) {
169 if ( m_references.at(refID).RefHasAlignments ) { break; }
172 m_currentRefID = refID;
174 m_isRegionSpecified = false;
176 return ( bam_seek(m_file, m_alignmentsBeginOffset, SEEK_SET) == 0 );
179 // get next alignment from specified region
180 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) {
182 // try to load 'next' read
183 if ( LoadNextAlignment(bAlignment) ) {
185 // if specified region, check for region overlap
186 if ( m_isRegionSpecified ) {
188 // if overlap, return true
189 if ( IsOverlap(bAlignment) ) { return true; }
190 // if not, try the next alignment
191 else { return GetNextAlignment(bAlignment); }
194 // not using region, valid read detected, return success
195 else { return true; }
198 // no valid alignment to load
202 void BamReader::SetCalculateAlignedBases(bool flag) {
203 m_isCalculateAlignedBases = flag;
206 int BamReader::BinsFromRegion(int refID, unsigned int left, uint16_t list[MAX_BIN]) {
208 // get region boundaries
209 uint32_t begin = left;
210 uint32_t end = m_references.at(refID).RefLength - 1;
212 // initialize list, bin '0' always a valid bin
216 // get rest of bins that contain this region
218 for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { list[i++] = k; }
219 for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { list[i++] = k; }
220 for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { list[i++] = k; }
221 for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { list[i++] = k; }
222 for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { list[i++] = k; }
224 // return number of bins stored
228 uint32_t BamReader::CalculateAlignmentEnd(const unsigned int& position, const vector<CigarOp>& cigarData) {
230 // initialize alignment end to starting position
231 uint32_t alignEnd = position;
233 // iterate over cigar operations
234 vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
235 vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
236 for ( ; cigarIter != cigarEnd; ++cigarIter) {
237 if ( (*cigarIter).Type == 'M' || (*cigarIter).Type == 'D' || (*cigarIter).Type == 'N') {
238 alignEnd += (*cigarIter).Length;
244 int64_t BamReader::GetOffset(int refID, unsigned int left) {
246 // make space for bins
247 uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);
249 // returns number of bins overlapping (left, right)
250 // stores indices of those bins in 'bins'
251 int numBins = BinsFromRegion(refID, left, bins);
253 // access index data for refID
254 RefIndex* refIndex = m_index->at(refID);
256 // get list of BAM bins for this reference sequence
257 BinVector* refBins = refIndex->first;
259 sort( refBins->begin(), refBins->end(), LookupKeyCompare<uint32_t, ChunkVector*>() );
261 // declare ChunkVector
262 ChunkVector regionChunks;
264 // declaure LinearOffsetVector
265 LinearOffsetVector* linearOffsets = refIndex->second;
267 // some sort of linear offset vs bin index voodoo... not completely sure what's going here
268 uint64_t minOffset = ((left>>BAM_LIDX_SHIFT) >= linearOffsets->size()) ? 0 : linearOffsets->at(left>>BAM_LIDX_SHIFT);
270 BinVector::iterator binBegin = refBins->begin();
271 BinVector::iterator binEnd = refBins->end();
273 // iterate over bins overlapping region, count chunks
274 for (int i = 0; i < numBins; ++i) {
276 // look for bin with ID=bin[i]
277 BinVector::iterator binIter = binBegin;
279 for ( ; binIter != binEnd; ++binIter ) {
281 // if found, increment n_off by number of chunks for each bin
282 if ( (*binIter).first == (uint32_t)bins[i] ) {
284 // iterate over chunks in that bin
285 ChunkVector* binChunks = (*binIter).second;
287 ChunkVector::iterator chunkIter = binChunks->begin();
288 ChunkVector::iterator chunkEnd = binChunks->end();
289 for ( ; chunkIter != chunkEnd; ++chunkIter) {
291 // if right bound of pair is greater than min_off (linear offset value), store pair
292 if ( (*chunkIter).second > minOffset) {
293 regionChunks.push_back( (*chunkIter) );
300 // clean up temp array
303 // there should be at least 1
304 if(regionChunks.size() > 0) { return -1; }
306 // sort chunks by start position
307 sort ( regionChunks.begin(), regionChunks.end(), LookupKeyCompare<uint64_t, uint64_t>() );
309 // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
310 int numOffsets = regionChunks.size();
311 for (int i = 1; i < numOffsets; ++i) {
312 if ( regionChunks.at(i-1).second >= regionChunks.at(i).first ) {
313 regionChunks.at(i-1).second = regionChunks.at(i).first;
317 // merge adjacent chunks
319 for (int i = 1; i < numOffsets; ++i) {
320 // if adjacent, expand boundaries of (merged) chunk
321 if ( (regionChunks.at(l).second>>16) == (regionChunks.at(i).first>>16) ) {
322 regionChunks.at(l).second = regionChunks.at(i).second;
324 // else, move on to next (merged) chunk index
325 else { regionChunks.at(++l) = regionChunks.at(i); }
328 // return beginning file offset of first chunk for region
329 return (int64_t)regionChunks.at(0).first;
332 bool BamReader::IsOverlap(BamAlignment& bAlignment) {
334 // if on different reference sequence, quit
335 if ( bAlignment.RefID != (unsigned int)m_currentRefID ) { return false; }
337 // read starts after left boundary
338 if ( bAlignment.Position >= m_currentLeft) { return true; }
341 uint32_t alignEnd = CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData);
343 // return whether alignment end overlaps left boundary
344 return ( alignEnd >= m_currentLeft );
347 bool BamReader::LoadHeader(void) {
349 // check to see if proper BAM header
351 if (bam_read(m_file, buf, 4) != 4) { return false; }
352 if (strncmp(buf, "BAM\001", 4)) {
353 fprintf(stderr, "wrong header type!\n");
357 // get BAM header text length
358 int32_t headerTextLength;
359 bam_read(m_file, &headerTextLength, 4);
361 // get BAM header text
362 char* headerText = (char*)calloc(headerTextLength + 1, 1);
363 bam_read(m_file, headerText, headerTextLength);
364 m_headerText = (string)((const char*)headerText);
366 // clean up calloc-ed temp variable
369 // get number of reference sequences
370 int32_t numberRefSeqs;
371 bam_read(m_file, &numberRefSeqs, 4);
372 if (numberRefSeqs == 0) { return false; }
374 m_references.reserve((int)numberRefSeqs);
376 // reference variables
377 int32_t refNameLength;
381 // iterate over all references in header
382 for (int i = 0; i != numberRefSeqs; ++i) {
384 // get length of reference name
385 bam_read(m_file, &refNameLength, 4);
386 refName = (char*)calloc(refNameLength, 1);
388 // get reference name and reference sequence length
389 bam_read(m_file, refName, refNameLength);
390 bam_read(m_file, &refLength, 4);
392 // store data for reference
394 aReference.RefName = (string)((const char*)refName);
395 aReference.RefLength = refLength;
396 m_references.push_back(aReference);
398 // clean up calloc-ed temp variable
405 bool BamReader::LoadIndex(void) {
407 // check to see if index file exists
409 if ( ( indexFile = fopen(m_indexFilename, "r") ) == 0 ) {
410 fprintf(stderr, "The alignment is not indexed. Please run SAMtools \'index\' command first.\n");
414 // see if index is valid BAM index
416 fread(magic, 1, 4, indexFile);
417 if (strncmp(magic, "BAI\1", 4)) {
418 fprintf(stderr, "Problem with index - wrong \'magic\' number.\n");
423 // get number of reference sequences
425 fread(&numRefSeqs, 4, 1, indexFile);
427 // intialize BamIndex data structure
428 m_index = new BamIndex;
429 m_index->reserve(numRefSeqs);
431 // iterate over reference sequences
432 for (unsigned int i = 0; i < numRefSeqs; ++i) {
434 // get number of bins for this reference sequence
436 fread(&numBins, 4, 1, indexFile);
438 if (numBins > 0) { m_references.at(i).RefHasAlignments = true; }
440 // intialize BinVector
441 BinVector* bins = new BinVector;
442 bins->reserve(numBins);
444 // iterate over bins for that reference sequence
445 for (int j = 0; j < numBins; ++j) {
449 fread(&binID, 4, 1, indexFile);
451 // get number of regionChunks in this bin
453 fread(&numChunks, 4, 1, indexFile);
455 // intialize ChunkVector
456 ChunkVector* regionChunks = new ChunkVector;
457 regionChunks->reserve(numChunks);
459 // iterate over regionChunks in this bin
460 for (unsigned int k = 0; k < numChunks; ++k) {
462 // get chunk boundaries (left, right)
465 fread(&left, 8, 1, indexFile);
466 fread(&right, 8, 1, indexFile);
469 regionChunks->push_back( ChunkPair(left, right) );
472 // save binID, chunkVector for this bin
473 bins->push_back( BamBin(binID, regionChunks) );
476 // load linear index for this reference sequence
478 // get number of linear offsets
479 int32_t numLinearOffsets;
480 fread(&numLinearOffsets, 4, 1, indexFile);
482 // intialize LinearOffsetVector
483 LinearOffsetVector* linearOffsets = new LinearOffsetVector;
484 linearOffsets->reserve(numLinearOffsets);
486 // iterate over linear offsets for this reference sequeence
487 for (int j = 0; j < numLinearOffsets; ++j) {
488 // get a linear offset
489 uint64_t linearOffset;
490 fread(&linearOffset, 8, 1, indexFile);
491 // store linear offset
492 linearOffsets->push_back(linearOffset);
495 // store index data for that reference sequence
496 m_index->push_back( new RefIndex(bins, linearOffsets) );
499 // close index file (.bai) and return
504 bool BamReader::LoadNextAlignment(BamAlignment& bAlignment) {
506 // check valid alignment block header data
511 int32_t bytesRead = 0;
513 // read in the 'block length' value, make sure it's not zero
514 if ( (ret = bam_read(m_file, &block_len, 4)) == 0 ) { return false; }
517 // read in core alignment data, make the right size of data was read
518 if ( bam_read(m_file, x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }
519 bytesRead += BAM_CORE_SIZE;
521 // set BamAlignment 'core' data
522 bAlignment.RefID = x[0];
523 bAlignment.Position = x[1];
524 bAlignment.Bin = x[2]>>16;
525 bAlignment.MapQuality = x[2]>>8&0xff;
526 bAlignment.AlignmentFlag = x[3]>>16;
527 bAlignment.MateRefID = x[5];
528 bAlignment.MatePosition = x[6];
529 bAlignment.InsertSize = x[7];
531 // fetch & store often-used lengths for character data parsing
532 unsigned int queryNameLength = x[2]&0xff;
533 unsigned int numCigarOperations = x[3]&0xffff;
534 unsigned int querySequenceLength = x[4];
536 // get length of character data
537 int dataLength = block_len - BAM_CORE_SIZE;
539 // set up destination buffers for character data
540 uint8_t* allCharData = (uint8_t*)calloc(sizeof(uint8_t), dataLength);
541 uint32_t* cigarData = (uint32_t*)(allCharData+queryNameLength);
543 const unsigned int tagDataOffset = (numCigarOperations * 4) + queryNameLength + (querySequenceLength + 1) / 2 + querySequenceLength;
544 const unsigned int tagDataLen = dataLength - tagDataOffset;
545 char* tagData = ((char*)allCharData) + tagDataOffset;
547 // get character data - make sure proper data size was read
548 if (bam_read(m_file, allCharData, dataLength) != dataLength) { return false; }
551 bytesRead += dataLength;
553 // clear out bases, qualities, aligned bases, CIGAR, and tag data
554 bAlignment.QueryBases.clear();
555 bAlignment.Qualities.clear();
556 bAlignment.AlignedBases.clear();
557 bAlignment.CigarData.clear();
558 bAlignment.TagData.clear();
561 bAlignment.Name = (string)((const char*)(allCharData));
565 uint8_t* s = ( allCharData + (numCigarOperations*4) + queryNameLength);
566 for (unsigned int i = 0; i < querySequenceLength; ++i) {
567 // numeric to char conversion
568 sprintf( singleBase, "%c", DNA_LOOKUP[ ( ( s[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ] );
569 // append character data to Bases
570 bAlignment.QueryBases.append( (const char*)singleBase );
573 // save sequence length
574 bAlignment.Length = bAlignment.QueryBases.length();
577 char singleQuality[4];
578 uint8_t* t = ( allCharData + (numCigarOperations*4) + queryNameLength + (querySequenceLength + 1)/2);
579 for (unsigned int i = 0; i < querySequenceLength; ++i) {
580 // numeric to char conversion
581 sprintf( singleQuality, "%c", ( t[i]+33 ) );
582 // append character data to Qualities
583 bAlignment.Qualities.append( (const char*)singleQuality );
586 // save CIGAR-related data;
588 for (unsigned int i = 0; i < numCigarOperations; ++i) {
590 // build CigarOp struct
592 op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
593 op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
596 bAlignment.CigarData.push_back(op);
598 // can skip this step if user wants to ignore
599 if (m_isCalculateAlignedBases) {
601 // build AlignedBases string
605 case ('I') : bAlignment.AlignedBases.append( bAlignment.QueryBases.substr(k, op.Length) ); // for 'M', 'I' - write bases
606 case ('S') : k += op.Length; // for 'S' - skip over query bases
610 case ('P') : bAlignment.AlignedBases.append( op.Length, '*' ); // for 'D', 'P' - write padding;
613 case ('N') : bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in query sequence
617 case ('H') : break; // for 'H' - do nothing, move to next op
619 default : assert(false); // shouldn't get here
625 // read in the tag data
626 bAlignment.TagData.resize(tagDataLen);
627 memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLen);
632 // (optional) read tag parsing data
637 // still data to read (auxiliary tags)
638 while ( bytesRead < block_len ) {
640 if ( bam_read(m_file, &data, 1) == 1 ) {
644 if (bytesRead == block_len && data != '\0') {
645 fprintf(stderr, "ERROR: Invalid termination of tag info at end of alignment block.\n");
650 if ( data == '\0' ) {
651 bAlignment.Tags.push_back(tag);
655 if ( (i == 1) && (i == 2) ) { tag.append(1, ':'); }
659 fprintf(stderr, "ERROR: Could not read tag info.\n");