-// BamReader.cpp
-
-// Derek Barnett
-// Marth Lab, Boston College
-// Last modified: 23 April 2009
-
-#include "BamReader.h"
-#include <iostream>
-using std::cerr;
-using std::endl;
-
-// static character constants
-const char* BamReader::DNA_LOOKUP = "=ACMGRSVTWYHKDBN";
-const char* BamReader::CIGAR_LOOKUP = "MIDNSHP";
-
-BamReader::BamReader(const char* filename, const char* indexFilename)
- : m_filename( (char*)filename )
- , m_indexFilename( (char*)indexFilename )
- , m_file(0)
- , m_index(NULL)
- , m_headerText("")
- , m_isOpen(false)
- , m_isIndexLoaded(false)
- , m_isRegionSpecified(false)
- , m_isCalculateAlignedBases(true)
- , m_currentRefID(0)
- , m_currentLeft(0)
- , m_alignmentsBeginOffset(0)
-{
- Open();
-}
-
-BamReader::~BamReader(void) {
-
- // close file
- if ( m_isOpen ) { Close(); }
-}
-
-// open BAM file
-bool BamReader::Open(void) {
-
- if (!m_isOpen && m_filename != NULL ) {
-
- // open file
- m_file = bam_open(m_filename, "r");
-
- // get header info && index info
- if ( (m_file != NULL) && LoadHeader() ) {
-
- // save file offset where alignments start
- m_alignmentsBeginOffset = bam_tell(m_file);
-
- // set open flag
- m_isOpen = true;
- }
-
- // try to open (and load) index data, if index file given
- if ( m_indexFilename != NULL ) {
- OpenIndex();
- }
- }
-
- return m_isOpen;
-}
-
-bool BamReader::OpenIndex(void) {
-
- if ( m_indexFilename && !m_isIndexLoaded ) {
- m_isIndexLoaded = LoadIndex();
- }
- return m_isIndexLoaded;
-}
-
-// close BAM file
-bool BamReader::Close(void) {
-
- if (m_isOpen) {
-
- // close file
- int ret = bam_close(m_file);
-
- // delete index info
- if ( m_index != NULL) { delete m_index; }
-
- // clear open flag
- m_isOpen = false;
-
- // clear index flag
- m_isIndexLoaded = false;
-
- // clear region flag
- m_isRegionSpecified = false;
-
- // return success/fail of bam_close
- return (ret == 0);
- }
-
- return true;
-}
-
-// get BAM filename
-const char* BamReader::Filename(void) const {
- return (const char*)m_filename;
-}
-
-// set BAM filename
-void BamReader::SetFilename(const char* filename) {
- m_filename = (char*)filename;
-}
-
-// get BAM Index filename
-const char* BamReader::IndexFilename(void) const {
- return (const char*)m_indexFilename;
-}
-
-// set BAM Index filename
-void BamReader::SetIndexFilename(const char* indexFilename) {
- m_indexFilename = (char*)indexFilename;
-}
-
-// return full header text
-const string BamReader::GetHeaderText(void) const {
- return m_headerText;
-}
-
-// return number of reference sequences in BAM file
-const int BamReader::GetReferenceCount(void) const {
- return m_references.size();
-}
-
-// get RefID from reference name
-const int BamReader::GetRefID(string refName) const {
-
- vector<string> refNames;
- RefVector::const_iterator refIter = m_references.begin();
- RefVector::const_iterator refEnd = m_references.end();
- for ( ; refIter != refEnd; ++refIter) {
- refNames.push_back( (*refIter).RefName );
- }
-
- // return 'index-of' refName (if not found, returns refNames.size())
- return Index( refNames.begin(), refNames.end(), refName );
-}
-
-const RefVector BamReader::GetReferenceData(void) const {
- return m_references;
-}
-
-bool BamReader::Jump(int refID, unsigned int left) {
-
- // if index available, and region is valid
- if ( m_isIndexLoaded && m_references.at(refID).RefHasAlignments && (left <= m_references.at(refID).RefLength) ) {
- m_currentRefID = refID;
- m_currentLeft = left;
- m_isRegionSpecified = true;
-
- int64_t offset = GetOffset(m_currentRefID, m_currentLeft);
- if ( offset == -1 ) { return false; }
- else { return ( bam_seek(m_file, offset, SEEK_SET) == 0 ); }
- }
- return false;
-}
-
-bool BamReader::Rewind(void) {
-
- int refID = 0;
- int refCount = m_references.size();
- for ( ; refID < refCount; ++refID ) {
- if ( m_references.at(refID).RefHasAlignments ) { break; }
- }
-
- m_currentRefID = refID;
- m_currentLeft = 0;
- m_isRegionSpecified = false;
-
- return ( bam_seek(m_file, m_alignmentsBeginOffset, SEEK_SET) == 0 );
-}
-
-// get next alignment from specified region
-bool BamReader::GetNextAlignment(BamAlignment& bAlignment) {
-
- // try to load 'next' read
- if ( LoadNextAlignment(bAlignment) ) {
-
- // if specified region, check for region overlap
- if ( m_isRegionSpecified ) {
-
- // if overlap, return true
- if ( IsOverlap(bAlignment) ) { return true; }
- // if not, try the next alignment
- else { return GetNextAlignment(bAlignment); }
- }
-
- // not using region, valid read detected, return success
- else { return true; }
- }
-
- // no valid alignment to load
- return false;
-}
-
-void BamReader::SetCalculateAlignedBases(bool flag) {
- m_isCalculateAlignedBases = flag;
-}
-
-int BamReader::BinsFromRegion(int refID, unsigned int left, uint16_t list[MAX_BIN]) {
-
- // get region boundaries
- uint32_t begin = left;
- uint32_t end = m_references.at(refID).RefLength - 1;
-
- // initialize list, bin '0' always a valid bin
- int i = 0;
- list[i++] = 0;
-
- // get rest of bins that contain this region
- unsigned int k;
- for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { list[i++] = k; }
- for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { list[i++] = k; }
- for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { list[i++] = k; }
- for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { list[i++] = k; }
- for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { list[i++] = k; }
-
- // return number of bins stored
- return i;
-}
-
-uint32_t BamReader::CalculateAlignmentEnd(const unsigned int& position, const vector<CigarOp>& cigarData) {
-
- // initialize alignment end to starting position
- uint32_t alignEnd = position;
-
- // iterate over cigar operations
- vector<CigarOp>::const_iterator cigarIter = cigarData.begin();
- vector<CigarOp>::const_iterator cigarEnd = cigarData.end();
- for ( ; cigarIter != cigarEnd; ++cigarIter) {
- if ( (*cigarIter).Type == 'M' || (*cigarIter).Type == 'D' || (*cigarIter).Type == 'N') {
- alignEnd += (*cigarIter).Length;
- }
- }
- return alignEnd;
-}
-
-int64_t BamReader::GetOffset(int refID, unsigned int left) {
-
- // make space for bins
- uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);
-
- // returns number of bins overlapping (left, right)
- // stores indices of those bins in 'bins'
- int numBins = BinsFromRegion(refID, left, bins);
-
- // access index data for refID
- RefIndex* refIndex = m_index->at(refID);
-
- // get list of BAM bins for this reference sequence
- BinVector* refBins = refIndex->first;
-
- sort( refBins->begin(), refBins->end(), LookupKeyCompare<uint32_t, ChunkVector*>() );
-
- // declare ChunkVector
- ChunkVector regionChunks;
-
- // declaure LinearOffsetVector
- LinearOffsetVector* linearOffsets = refIndex->second;
-
- // some sort of linear offset vs bin index voodoo... not completely sure what's going here
- uint64_t minOffset = ((left>>BAM_LIDX_SHIFT) >= linearOffsets->size()) ? 0 : linearOffsets->at(left>>BAM_LIDX_SHIFT);
-
- BinVector::iterator binBegin = refBins->begin();
- BinVector::iterator binEnd = refBins->end();
-
- // iterate over bins overlapping region, count chunks
- for (int i = 0; i < numBins; ++i) {
-
- // look for bin with ID=bin[i]
- BinVector::iterator binIter = binBegin;
-
- for ( ; binIter != binEnd; ++binIter ) {
-
- // if found, increment n_off by number of chunks for each bin
- if ( (*binIter).first == (uint32_t)bins[i] ) {
-
- // iterate over chunks in that bin
- ChunkVector* binChunks = (*binIter).second;
-
- ChunkVector::iterator chunkIter = binChunks->begin();
- ChunkVector::iterator chunkEnd = binChunks->end();
- for ( ; chunkIter != chunkEnd; ++chunkIter) {
-
- // if right bound of pair is greater than min_off (linear offset value), store pair
- if ( (*chunkIter).second > minOffset) {
- regionChunks.push_back( (*chunkIter) );
- }
- }
- }
- }
- }
-
- // clean up temp array
- free(bins);
-
- // there should be at least 1
- if(regionChunks.size() > 0) { return -1; }
-
- // sort chunks by start position
- sort ( regionChunks.begin(), regionChunks.end(), LookupKeyCompare<uint64_t, uint64_t>() );
-
- // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
- int numOffsets = regionChunks.size();
- for (int i = 1; i < numOffsets; ++i) {
- if ( regionChunks.at(i-1).second >= regionChunks.at(i).first ) {
- regionChunks.at(i-1).second = regionChunks.at(i).first;
- }
- }
-
- // merge adjacent chunks
- int l = 0;
- for (int i = 1; i < numOffsets; ++i) {
- // if adjacent, expand boundaries of (merged) chunk
- if ( (regionChunks.at(l).second>>16) == (regionChunks.at(i).first>>16) ) {
- regionChunks.at(l).second = regionChunks.at(i).second;
- }
- // else, move on to next (merged) chunk index
- else { regionChunks.at(++l) = regionChunks.at(i); }
- }
-
- // return beginning file offset of first chunk for region
- return (int64_t)regionChunks.at(0).first;
-}
-
-bool BamReader::IsOverlap(BamAlignment& bAlignment) {
-
- // if on different reference sequence, quit
- if ( bAlignment.RefID != (unsigned int)m_currentRefID ) { return false; }
-
- // read starts after left boundary
- if ( bAlignment.Position >= m_currentLeft) { return true; }
-
- // get alignment end
- uint32_t alignEnd = CalculateAlignmentEnd(bAlignment.Position, bAlignment.CigarData);
-
- // return whether alignment end overlaps left boundary
- return ( alignEnd >= m_currentLeft );
-}
-
-bool BamReader::LoadHeader(void) {
-
- // check to see if proper BAM header
- char buf[4];
- if (bam_read(m_file, buf, 4) != 4) { return false; }
- if (strncmp(buf, "BAM\001", 4)) {
- fprintf(stderr, "wrong header type!\n");
- return false;
- }
-
- // get BAM header text length
- int32_t headerTextLength;
- bam_read(m_file, &headerTextLength, 4);
-
- // get BAM header text
- char* headerText = (char*)calloc(headerTextLength + 1, 1);
- bam_read(m_file, headerText, headerTextLength);
- m_headerText = (string)((const char*)headerText);
-
- // clean up calloc-ed temp variable
- free(headerText);
-
- // get number of reference sequences
- int32_t numberRefSeqs;
- bam_read(m_file, &numberRefSeqs, 4);
- if (numberRefSeqs == 0) { return false; }
-
- m_references.reserve((int)numberRefSeqs);
-
- // reference variables
- int32_t refNameLength;
- char* refName;
- uint32_t refLength;
-
- // iterate over all references in header
- for (int i = 0; i != numberRefSeqs; ++i) {
-
- // get length of reference name
- bam_read(m_file, &refNameLength, 4);
- refName = (char*)calloc(refNameLength, 1);
-
- // get reference name and reference sequence length
- bam_read(m_file, refName, refNameLength);
- bam_read(m_file, &refLength, 4);
-
- // store data for reference
- RefData aReference;
- aReference.RefName = (string)((const char*)refName);
- aReference.RefLength = refLength;
- m_references.push_back(aReference);
-
- // clean up calloc-ed temp variable
- free(refName);
- }
-
- return true;
-}
-
-bool BamReader::LoadIndex(void) {
-
- // check to see if index file exists
- FILE* indexFile;
- if ( ( indexFile = fopen(m_indexFilename, "r") ) == 0 ) {
- fprintf(stderr, "The alignment is not indexed. Please run SAMtools \'index\' command first.\n");
- return false;
- }
-
- // see if index is valid BAM index
- char magic[4];
- fread(magic, 1, 4, indexFile);
- if (strncmp(magic, "BAI\1", 4)) {
- fprintf(stderr, "Problem with index - wrong \'magic\' number.\n");
- fclose(indexFile);
- return false;
- }
-
- // get number of reference sequences
- uint32_t numRefSeqs;
- fread(&numRefSeqs, 4, 1, indexFile);
-
- // intialize BamIndex data structure
- m_index = new BamIndex;
- m_index->reserve(numRefSeqs);
-
- // iterate over reference sequences
- for (unsigned int i = 0; i < numRefSeqs; ++i) {
-
- // get number of bins for this reference sequence
- int32_t numBins;
- fread(&numBins, 4, 1, indexFile);
-
- if (numBins > 0) { m_references.at(i).RefHasAlignments = true; }
-
- // intialize BinVector
- BinVector* bins = new BinVector;
- bins->reserve(numBins);
-
- // iterate over bins for that reference sequence
- for (int j = 0; j < numBins; ++j) {
-
- // get binID
- uint32_t binID;
- fread(&binID, 4, 1, indexFile);
-
- // get number of regionChunks in this bin
- uint32_t numChunks;
- fread(&numChunks, 4, 1, indexFile);
-
- // intialize ChunkVector
- ChunkVector* regionChunks = new ChunkVector;
- regionChunks->reserve(numChunks);
-
- // iterate over regionChunks in this bin
- for (unsigned int k = 0; k < numChunks; ++k) {
-
- // get chunk boundaries (left, right)
- uint64_t left;
- uint64_t right;
- fread(&left, 8, 1, indexFile);
- fread(&right, 8, 1, indexFile);
-
- // save ChunkPair
- regionChunks->push_back( ChunkPair(left, right) );
- }
-
- // save binID, chunkVector for this bin
- bins->push_back( BamBin(binID, regionChunks) );
- }
-
- // load linear index for this reference sequence
-
- // get number of linear offsets
- int32_t numLinearOffsets;
- fread(&numLinearOffsets, 4, 1, indexFile);
-
- // intialize LinearOffsetVector
- LinearOffsetVector* linearOffsets = new LinearOffsetVector;
- linearOffsets->reserve(numLinearOffsets);
-
- // iterate over linear offsets for this reference sequeence
- for (int j = 0; j < numLinearOffsets; ++j) {
- // get a linear offset
- uint64_t linearOffset;
- fread(&linearOffset, 8, 1, indexFile);
- // store linear offset
- linearOffsets->push_back(linearOffset);
- }
-
- // store index data for that reference sequence
- m_index->push_back( new RefIndex(bins, linearOffsets) );
- }
-
- // close index file (.bai) and return
- fclose(indexFile);
- return true;
-}
-
-bool BamReader::LoadNextAlignment(BamAlignment& bAlignment) {
-
- // check valid alignment block header data
- int32_t block_len;
- int32_t ret;
- uint32_t x[8];
-
- int32_t bytesRead = 0;
-
- // read in the 'block length' value, make sure it's not zero
- if ( (ret = bam_read(m_file, &block_len, 4)) == 0 ) { return false; }
- bytesRead += 4;
-
- // read in core alignment data, make the right size of data was read
- if ( bam_read(m_file, x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }
- bytesRead += BAM_CORE_SIZE;
-
- // set BamAlignment 'core' data
- bAlignment.RefID = x[0];
- bAlignment.Position = x[1];
- bAlignment.Bin = x[2]>>16;
- bAlignment.MapQuality = x[2]>>8&0xff;
- bAlignment.AlignmentFlag = x[3]>>16;
- bAlignment.MateRefID = x[5];
- bAlignment.MatePosition = x[6];
- bAlignment.InsertSize = x[7];
-
- // fetch & store often-used lengths for character data parsing
- unsigned int queryNameLength = x[2]&0xff;
- unsigned int numCigarOperations = x[3]&0xffff;
- unsigned int querySequenceLength = x[4];
-
- // get length of character data
- int dataLength = block_len - BAM_CORE_SIZE;
-
- // set up destination buffers for character data
- uint8_t* allCharData = (uint8_t*)calloc(sizeof(uint8_t), dataLength);
- uint32_t* cigarData = (uint32_t*)(allCharData+queryNameLength);
-
- const unsigned int tagDataOffset = (numCigarOperations * 4) + queryNameLength + (querySequenceLength + 1) / 2 + querySequenceLength;
- const unsigned int tagDataLen = dataLength - tagDataOffset;
- char* tagData = ((char*)allCharData) + tagDataOffset;
-
- // get character data - make sure proper data size was read
- if (bam_read(m_file, allCharData, dataLength) != dataLength) { return false; }
- else {
-
- bytesRead += dataLength;
-
- // clear out bases, qualities, aligned bases, CIGAR, and tag data
- bAlignment.QueryBases.clear();
- bAlignment.Qualities.clear();
- bAlignment.AlignedBases.clear();
- bAlignment.CigarData.clear();
- bAlignment.TagData.clear();
-
- // save name
- bAlignment.Name = (string)((const char*)(allCharData));
-
- // save bases
- char singleBase[2];
- uint8_t* s = ( allCharData + (numCigarOperations*4) + queryNameLength);
- for (unsigned int i = 0; i < querySequenceLength; ++i) {
- // numeric to char conversion
- sprintf( singleBase, "%c", DNA_LOOKUP[ ( ( s[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ] );
- // append character data to Bases
- bAlignment.QueryBases.append( (const char*)singleBase );
- }
-
- // save sequence length
- bAlignment.Length = bAlignment.QueryBases.length();
-
- // save qualities
- char singleQuality[4];
- uint8_t* t = ( allCharData + (numCigarOperations*4) + queryNameLength + (querySequenceLength + 1)/2);
- for (unsigned int i = 0; i < querySequenceLength; ++i) {
- // numeric to char conversion
- sprintf( singleQuality, "%c", ( t[i]+33 ) );
- // append character data to Qualities
- bAlignment.Qualities.append( (const char*)singleQuality );
- }
-
- // save CIGAR-related data;
- int k = 0;
- for (unsigned int i = 0; i < numCigarOperations; ++i) {
-
- // build CigarOp struct
- CigarOp op;
- op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
- op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
-
- // save CigarOp
- bAlignment.CigarData.push_back(op);
-
- // can skip this step if user wants to ignore
- if (m_isCalculateAlignedBases) {
-
- // build AlignedBases string
- switch (op.Type) {
-
- case ('M') :
- case ('I') : bAlignment.AlignedBases.append( bAlignment.QueryBases.substr(k, op.Length) ); // for 'M', 'I' - write bases
- case ('S') : k += op.Length; // for 'S' - skip over query bases
- break;
-
- case ('D') :
- case ('P') : bAlignment.AlignedBases.append( op.Length, '*' ); // for 'D', 'P' - write padding;
- break;
-
- case ('N') : bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in query sequence
- k += op.Length;
- break;
-
- case ('H') : break; // for 'H' - do nothing, move to next op
-
- default : assert(false); // shouldn't get here
- break;
- }
- }
- }
-
- // read in the tag data
- bAlignment.TagData.resize(tagDataLen);
- memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLen);
- }
- free(allCharData);
-
-/*
- // (optional) read tag parsing data
- string tag;
- char data;
- int i = 0;
-
- // still data to read (auxiliary tags)
- while ( bytesRead < block_len ) {
-
- if ( bam_read(m_file, &data, 1) == 1 ) {
-
- ++bytesRead;
-
- if (bytesRead == block_len && data != '\0') {
- fprintf(stderr, "ERROR: Invalid termination of tag info at end of alignment block.\n");
- exit(1);
- }
-
- tag.append(1, data);
- if ( data == '\0' ) {
- bAlignment.Tags.push_back(tag);
- tag = "";
- i = 0;
- } else {
- if ( (i == 1) && (i == 2) ) { tag.append(1, ':'); }
- ++i;
- }
- } else {
- fprintf(stderr, "ERROR: Could not read tag info.\n");
- exit(1);
- }
- }
-*/
- return true;
-}
+// ***************************************************************************\r
+// BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg\r
+// Marth Lab, Department of Biology, Boston College\r
+// All rights reserved.\r
+// ---------------------------------------------------------------------------\r
+// Last modified: 15 July 2010 (DB)\r
+// ---------------------------------------------------------------------------\r
+// Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
+// Institute.\r
+// ---------------------------------------------------------------------------\r
+// Provides the basic functionality for reading BAM files\r
+// ***************************************************************************\r
+\r
+// C++ includes\r
+#include <algorithm>\r
+#include <iterator>\r
+#include <string>\r
+#include <vector>\r
+#include <iostream>\r
+\r
+// BamTools includes\r
+#include "BGZF.h"\r
+#include "BamReader.h"\r
+#include "BamIndex.h"\r
+using namespace BamTools;\r
+using namespace std;\r
+\r
+struct BamReader::BamReaderPrivate {\r
+\r
+ // -------------------------------\r
+ // structs, enums, typedefs\r
+ // -------------------------------\r
+ enum RegionState { BEFORE_REGION = 0\r
+ , WITHIN_REGION\r
+ , AFTER_REGION\r
+ };\r
+\r
+ // -------------------------------\r
+ // data members\r
+ // -------------------------------\r
+\r
+ // general file data\r
+ BgzfData mBGZF;\r
+ string HeaderText;\r
+ //BamIndex Index;\r
+ BamIndex* NewIndex;\r
+ RefVector References;\r
+ bool IsIndexLoaded;\r
+ int64_t AlignmentsBeginOffset;\r
+ string Filename;\r
+ string IndexFilename;\r
+ \r
+ // system data\r
+ bool IsBigEndian;\r
+\r
+ // user-specified region values\r
+ BamRegion Region;\r
+ bool IsLeftBoundSpecified;\r
+ bool IsRightBoundSpecified;\r
+ \r
+ bool IsRegionSpecified;\r
+ int CurrentRefID;\r
+ int CurrentLeft;\r
+\r
+ // parent BamReader\r
+ BamReader* Parent;\r
+ \r
+ // BAM character constants\r
+ const char* DNA_LOOKUP;\r
+ const char* CIGAR_LOOKUP;\r
+\r
+ // -------------------------------\r
+ // constructor & destructor\r
+ // -------------------------------\r
+ BamReaderPrivate(BamReader* parent);\r
+ ~BamReaderPrivate(void);\r
+\r
+ // -------------------------------\r
+ // "public" interface\r
+ // -------------------------------\r
+\r
+ // file operations\r
+ void Close(void);\r
+ bool Jump(int refID, int position = 0);\r
+ bool Open(const string& filename, const string& indexFilename = "");\r
+ bool Rewind(void);\r
+ bool SetRegion(const BamRegion& region);\r
+\r
+ // access alignment data\r
+ bool GetNextAlignment(BamAlignment& bAlignment);\r
+ bool GetNextAlignmentCore(BamAlignment& bAlignment);\r
+\r
+ // access auxiliary data\r
+ int GetReferenceID(const string& refName) const;\r
+\r
+ // index operations\r
+ bool CreateIndex(bool useDefaultIndex);\r
+\r
+ // -------------------------------\r
+ // internal methods\r
+ // -------------------------------\r
+\r
+ // *** reading alignments and auxiliary data *** //\r
+\r
+ // fills out character data for BamAlignment data\r
+ bool BuildCharData(BamAlignment& bAlignment);\r
+ // checks to see if alignment overlaps current region\r
+ RegionState IsOverlap(BamAlignment& bAlignment);\r
+ // retrieves header text from BAM file\r
+ void LoadHeaderData(void);\r
+ // retrieves BAM alignment under file pointer\r
+ bool LoadNextAlignment(BamAlignment& bAlignment);\r
+ // builds reference data structure from BAM file\r
+ void LoadReferenceData(void);\r
+\r
+ // *** index file handling *** //\r
+\r
+ // clear out inernal index data structure\r
+ void ClearIndex(void);\r
+ // loads index from BAM index file\r
+ bool LoadIndex(void);\r
+};\r
+\r
+// -----------------------------------------------------\r
+// BamReader implementation (wrapper around BRPrivate)\r
+// -----------------------------------------------------\r
+// constructor\r
+BamReader::BamReader(void) {\r
+ d = new BamReaderPrivate(this);\r
+}\r
+\r
+// destructor\r
+BamReader::~BamReader(void) {\r
+ delete d;\r
+ d = 0;\r
+}\r
+\r
+// file operations\r
+void BamReader::Close(void) { d->Close(); }\r
+bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }\r
+bool BamReader::Jump(int refID, int position) { \r
+ d->Region.LeftRefID = refID;\r
+ d->Region.LeftPosition = position;\r
+ d->IsLeftBoundSpecified = true;\r
+ d->IsRightBoundSpecified = false;\r
+ return d->Jump(refID, position); \r
+}\r
+bool BamReader::Open(const string& filename, const string& indexFilename) { return d->Open(filename, indexFilename); }\r
+bool BamReader::Rewind(void) { return d->Rewind(); }\r
+bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }\r
+bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {\r
+ return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );\r
+}\r
+\r
+// access alignment data\r
+bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }\r
+bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }\r
+\r
+// access auxiliary data\r
+const string BamReader::GetHeaderText(void) const { return d->HeaderText; }\r
+int BamReader::GetReferenceCount(void) const { return d->References.size(); }\r
+const RefVector& BamReader::GetReferenceData(void) const { return d->References; }\r
+int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }\r
+const std::string BamReader::GetFilename(void) const { return d->Filename; }\r
+\r
+// index operations\r
+bool BamReader::CreateIndex(bool useDefaultIndex) { return d->CreateIndex(useDefaultIndex); }\r
+\r
+// -----------------------------------------------------\r
+// BamReaderPrivate implementation\r
+// -----------------------------------------------------\r
+\r
+// constructor\r
+BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)\r
+ : NewIndex(0)\r
+ , IsIndexLoaded(false)\r
+ , AlignmentsBeginOffset(0)\r
+ , IsLeftBoundSpecified(false)\r
+ , IsRightBoundSpecified(false)\r
+ , IsRegionSpecified(false)\r
+ , CurrentRefID(0)\r
+ , CurrentLeft(0)\r
+ , Parent(parent)\r
+ , DNA_LOOKUP("=ACMGRSVTWYHKDBN")\r
+ , CIGAR_LOOKUP("MIDNSHP")\r
+{ \r
+ IsBigEndian = SystemIsBigEndian();\r
+}\r
+\r
+// destructor\r
+BamReader::BamReaderPrivate::~BamReaderPrivate(void) {\r
+ Close();\r
+}\r
+\r
+bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {\r
+ \r
+ // calculate character lengths/offsets\r
+ const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;\r
+ const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;\r
+ const unsigned int seqDataOffset = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);\r
+ const unsigned int qualDataOffset = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;\r
+ const unsigned int tagDataOffset = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;\r
+ const unsigned int tagDataLength = dataLength - tagDataOffset;\r
+ \r
+ // set up char buffers\r
+ const char* allCharData = bAlignment.SupportData.AllCharData.data();\r
+ uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);\r
+ const char* seqData = ((const char*)allCharData) + seqDataOffset;\r
+ const char* qualData = ((const char*)allCharData) + qualDataOffset;\r
+ char* tagData = ((char*)allCharData) + tagDataOffset;\r
+ \r
+ // store alignment name (depends on null char as terminator)\r
+ bAlignment.Name.assign((const char*)(allCharData)); \r
+ \r
+ // save CigarOps \r
+ CigarOp op;\r
+ bAlignment.CigarData.clear();\r
+ bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);\r
+ for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {\r
+\r
+ // swap if necessary\r
+ if ( IsBigEndian ) { SwapEndian_32(cigarData[i]); }\r
+ \r
+ // build CigarOp structure\r
+ op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);\r
+ op.Type = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];\r
+\r
+ // save CigarOp\r
+ bAlignment.CigarData.push_back(op);\r
+ }\r
+ \r
+ \r
+ // save query sequence\r
+ bAlignment.QueryBases.clear();\r
+ bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);\r
+ for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {\r
+ char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];\r
+ bAlignment.QueryBases.append(1, singleBase);\r
+ }\r
+ \r
+ // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character\r
+ bAlignment.Qualities.clear();\r
+ bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);\r
+ for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {\r
+ char singleQuality = (char)(qualData[i]+33);\r
+ bAlignment.Qualities.append(1, singleQuality);\r
+ }\r
+ \r
+ // if QueryBases is empty (and this is a allowed case)\r
+ if ( bAlignment.QueryBases.empty() ) \r
+ bAlignment.AlignedBases = bAlignment.QueryBases;\r
+ \r
+ // if QueryBases contains data, then build AlignedBases using CIGAR data\r
+ else {\r
+ \r
+ // resize AlignedBases\r
+ bAlignment.AlignedBases.clear();\r
+ bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);\r
+ \r
+ // iterate over CigarOps\r
+ int k = 0;\r
+ vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();\r
+ vector<CigarOp>::const_iterator cigarEnd = bAlignment.CigarData.end();\r
+ for ( ; cigarIter != cigarEnd; ++cigarIter ) {\r
+ \r
+ const CigarOp& op = (*cigarIter);\r
+ switch(op.Type) {\r
+ \r
+ case ('M') :\r
+ case ('I') :\r
+ bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases\r
+ // fall through\r
+ \r
+ case ('S') :\r
+ k += op.Length; // for 'S' - soft clip, skip over query bases\r
+ break;\r
+ \r
+ case ('D') :\r
+ bAlignment.AlignedBases.append(op.Length, '-'); // for 'D' - write gap character\r
+ break;\r
+ \r
+ case ('P') :\r
+ bAlignment.AlignedBases.append( op.Length, '*' ); // for 'P' - write padding character\r
+ break;\r
+ \r
+ case ('N') :\r
+ bAlignment.AlignedBases.append( op.Length, 'N' ); // for 'N' - write N's, skip bases in original query sequence\r
+ break;\r
+ \r
+ case ('H') :\r
+ break; // for 'H' - hard clip, do nothing to AlignedBases, move to next op\r
+ \r
+ default:\r
+ printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here\r
+ exit(1);\r
+ }\r
+ }\r
+ }\r
+ \r
+ // -----------------------\r
+ // Added: 3-25-2010 DB\r
+ // Fixed: endian-correctness for tag data\r
+ // -----------------------\r
+ if ( IsBigEndian ) {\r
+ int i = 0;\r
+ while ( (unsigned int)i < tagDataLength ) {\r
+ \r
+ i += 2; // skip tag type (e.g. "RG", "NM", etc)\r
+ uint8_t type = toupper(tagData[i]); // lower & upper case letters have same meaning \r
+ ++i; // skip value type\r
+ \r
+ switch (type) {\r
+ \r
+ case('A') :\r
+ case('C') : \r
+ ++i;\r
+ break;\r
+\r
+ case('S') : \r
+ SwapEndian_16p(&tagData[i]); \r
+ i += sizeof(uint16_t);\r
+ break;\r
+ \r
+ case('F') :\r
+ case('I') : \r
+ SwapEndian_32p(&tagData[i]);\r
+ i += sizeof(uint32_t);\r
+ break;\r
+ \r
+ case('D') : \r
+ SwapEndian_64p(&tagData[i]);\r
+ i += sizeof(uint64_t);\r
+ break;\r
+ \r
+ case('H') :\r
+ case('Z') : \r
+ while (tagData[i]) { ++i; }\r
+ ++i; // increment one more for null terminator\r
+ break;\r
+ \r
+ default : \r
+ printf("ERROR: Invalid tag value type\n"); // shouldn't get here\r
+ exit(1);\r
+ }\r
+ }\r
+ }\r
+ \r
+ // store TagData\r
+ bAlignment.TagData.clear();\r
+ bAlignment.TagData.resize(tagDataLength);\r
+ memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);\r
+ \r
+ // clear the core-only flag\r
+ bAlignment.SupportData.HasCoreOnly = false;\r
+ \r
+ // return success\r
+ return true;\r
+}\r
+\r
+// clear index data structure\r
+void BamReader::BamReaderPrivate::ClearIndex(void) {\r
+ delete NewIndex;\r
+ NewIndex = 0;\r
+}\r
+\r
+// closes the BAM file\r
+void BamReader::BamReaderPrivate::Close(void) {\r
+ \r
+ // close BGZF file stream\r
+ mBGZF.Close();\r
+ \r
+ // clear out index data\r
+ ClearIndex();\r
+ \r
+ // clear out header data\r
+ HeaderText.clear();\r
+ \r
+ // clear out region flags\r
+ IsLeftBoundSpecified = false;\r
+ IsRightBoundSpecified = false;\r
+ IsRegionSpecified = false;\r
+}\r
+\r
+// create BAM index from BAM file (keep structure in memory) and write to default index output file\r
+bool BamReader::BamReaderPrivate::CreateIndex(bool useDefaultIndex) {\r
+\r
+ // clear out prior index data\r
+ ClearIndex();\r
+ \r
+ // create default index\r
+ if ( useDefaultIndex )\r
+ NewIndex = new BamDefaultIndex(&mBGZF, Parent, IsBigEndian);\r
+ // create BamTools 'custom' index\r
+ else\r
+ NewIndex = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);\r
+ \r
+ bool ok = true;\r
+ ok &= NewIndex->Build();\r
+ ok &= NewIndex->Write(Filename); \r
+ \r
+ // return success/fail\r
+ return ok;\r
+}\r
+\r
+// get next alignment (from specified region, if given)\r
+bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {\r
+\r
+ // if valid alignment found, attempt to parse char data, and return success/failure\r
+ if ( GetNextAlignmentCore(bAlignment) )\r
+ return BuildCharData(bAlignment);\r
+ \r
+ // no valid alignment found\r
+ else\r
+ return false;\r
+}\r
+\r
+// retrieves next available alignment core data (returns success/fail)\r
+// ** DOES NOT parse any character data (read name, bases, qualities, tag data)\r
+// these can be accessed, if necessary, from the supportData \r
+// useful for operations requiring ONLY positional or other alignment-related information\r
+bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {\r
+\r
+ // if valid alignment available\r
+ if ( LoadNextAlignment(bAlignment) ) {\r
+\r
+ // set core-only flag\r
+ bAlignment.SupportData.HasCoreOnly = true;\r
+ \r
+ // if region not specified, return success\r
+ if ( !IsLeftBoundSpecified ) return true;\r
+\r
+ // determine region state (before, within, after)\r
+ BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);\r
+ \r
+ // if alignment lies after region, return false\r
+ if ( state == AFTER_REGION ) \r
+ return false;\r
+\r
+ while ( state != WITHIN_REGION ) {\r
+ // if no valid alignment available (likely EOF) return failure\r
+ if ( !LoadNextAlignment(bAlignment) ) return false;\r
+ // if alignment lies after region, return false (no available read within region)\r
+ state = IsOverlap(bAlignment);\r
+ if ( state == AFTER_REGION) return false;\r
+ \r
+ }\r
+\r
+ // return success (alignment found that overlaps region)\r
+ return true;\r
+ }\r
+\r
+ // no valid alignment\r
+ else\r
+ return false;\r
+}\r
+\r
+// returns RefID for given RefName (returns References.size() if not found)\r
+int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {\r
+\r
+ // retrieve names from reference data\r
+ vector<string> refNames;\r
+ RefVector::const_iterator refIter = References.begin();\r
+ RefVector::const_iterator refEnd = References.end();\r
+ for ( ; refIter != refEnd; ++refIter) {\r
+ refNames.push_back( (*refIter).RefName );\r
+ }\r
+\r
+ // return 'index-of' refName ( if not found, returns refNames.size() )\r
+ return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));\r
+}\r
+\r
+// returns region state - whether alignment ends before, overlaps, or starts after currently specified region\r
+// this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true\r
+BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {\r
+ \r
+ // --------------------------------------------------\r
+ // check alignment start against right bound cutoff\r
+ \r
+ // if full region of interest was given\r
+ if ( IsRightBoundSpecified ) {\r
+ \r
+ // read starts on right bound reference, but AFTER right bound position\r
+ if ( bAlignment.RefID == Region.RightRefID && bAlignment.Position > Region.RightPosition )\r
+ return AFTER_REGION;\r
+ \r
+ // if read starts on reference AFTER right bound, return false\r
+ if ( bAlignment.RefID > Region.RightRefID ) \r
+ return AFTER_REGION;\r
+ }\r
+ \r
+ // --------------------------------------------------------\r
+ // no right bound given OR read starts before right bound\r
+ // so, check if it overlaps left bound \r
+ \r
+ // if read starts on left bound reference AND after left boundary, return success\r
+ if ( bAlignment.RefID == Region.LeftRefID && bAlignment.Position >= Region.LeftPosition)\r
+ return WITHIN_REGION;\r
+ \r
+ // if read is on any reference sequence before left bound, return false\r
+ if ( bAlignment.RefID < Region.LeftRefID )\r
+ return BEFORE_REGION;\r
+\r
+ // --------------------------------------------------------\r
+ // read is on left bound reference, but starts before left bound position\r
+\r
+ // if it overlaps, return WITHIN_REGION\r
+ if ( bAlignment.GetEndPosition() >= Region.LeftPosition )\r
+ return WITHIN_REGION;\r
+ // else begins before left bound position\r
+ else\r
+ return BEFORE_REGION;\r
+}\r
+\r
+// jumps to specified region(refID, leftBound) in BAM file, returns success/fail\r
+bool BamReader::BamReaderPrivate::Jump(int refID, int position) {\r
+\r
+ // -----------------------------------------------------------------------\r
+ // check for existing index \r
+ if ( NewIndex == 0 ) return false; \r
+ // see if reference has alignments\r
+ if ( !NewIndex->HasAlignments(refID) ) return false; \r
+ // make sure position is valid\r
+ if ( position > References.at(refID).RefLength ) return false;\r
+ \r
+ // determine possible offsets\r
+ vector<int64_t> offsets;\r
+ if ( !NewIndex->GetOffsets(Region, IsRightBoundSpecified, offsets) ) {\r
+ printf("ERROR: Could not jump: unable to calculate offset for specified region.\n");\r
+ return false;\r
+ }\r
+ \r
+ // iterate through offsets\r
+ BamAlignment bAlignment;\r
+ bool result = true;\r
+ for ( vector<int64_t>::const_iterator o = offsets.begin(); o != offsets.end(); ++o) {\r
+ \r
+ // attempt seek & load first available alignment\r
+ result &= mBGZF.Seek(*o);\r
+ LoadNextAlignment(bAlignment);\r
+ \r
+ // if this alignment corresponds to desired position\r
+ // return success of seeking back to 'current offset'\r
+ if ( (bAlignment.RefID == refID && bAlignment.Position + bAlignment.Length > position) || (bAlignment.RefID > refID) ) {\r
+ if ( o != offsets.begin() ) --o;\r
+ return mBGZF.Seek(*o);\r
+ }\r
+ }\r
+ \r
+ return result;\r
+}\r
+\r
+// load BAM header data\r
+void BamReader::BamReaderPrivate::LoadHeaderData(void) {\r
+\r
+ // check to see if proper BAM header\r
+ char buffer[4];\r
+ if (mBGZF.Read(buffer, 4) != 4) {\r
+ printf("Could not read header type\n");\r
+ exit(1);\r
+ }\r
+\r
+ if (strncmp(buffer, "BAM\001", 4)) {\r
+ printf("wrong header type!\n");\r
+ exit(1);\r
+ }\r
+\r
+ // get BAM header text length\r
+ mBGZF.Read(buffer, 4);\r
+ unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);\r
+ if ( IsBigEndian ) { SwapEndian_32(headerTextLength); }\r
+ \r
+ // get BAM header text\r
+ char* headerText = (char*)calloc(headerTextLength + 1, 1);\r
+ mBGZF.Read(headerText, headerTextLength);\r
+ HeaderText = (string)((const char*)headerText);\r
+\r
+ // clean up calloc-ed temp variable\r
+ free(headerText);\r
+}\r
+\r
+// load existing index data from BAM index file (".bai"), return success/fail\r
+bool BamReader::BamReaderPrivate::LoadIndex(void) {\r
+\r
+ // clear out any existing index data\r
+ ClearIndex();\r
+\r
+ // skip if index file empty\r
+ if ( IndexFilename.empty() )\r
+ return false;\r
+\r
+ // check supplied filename for index type\r
+ size_t defaultExtensionFound = IndexFilename.find(".bai");\r
+ size_t customExtensionFound = IndexFilename.find(".bti");\r
+ \r
+ // if SAM/BAM default (".bai")\r
+ if ( defaultExtensionFound != string::npos )\r
+ NewIndex = new BamDefaultIndex(&mBGZF, Parent, IsBigEndian);\r
+ \r
+ // if BamTools custom index (".bti")\r
+ else if ( customExtensionFound != string::npos )\r
+ NewIndex = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);\r
+ \r
+ // else unknown\r
+ else {\r
+ printf("ERROR: Unknown index file extension.\n");\r
+ return false;\r
+ }\r
+ \r
+ // return success of loading index data\r
+ return NewIndex->Load(IndexFilename);\r
+}\r
+\r
+// populates BamAlignment with alignment data under file pointer, returns success/fail\r
+bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {\r
+\r
+ // read in the 'block length' value, make sure it's not zero\r
+ char buffer[4];\r
+ mBGZF.Read(buffer, 4);\r
+ bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);\r
+ if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }\r
+ if ( bAlignment.SupportData.BlockLength == 0 ) { return false; }\r
+\r
+ // read in core alignment data, make sure the right size of data was read\r
+ char x[BAM_CORE_SIZE];\r
+ if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) { return false; }\r
+\r
+ if ( IsBigEndian ) {\r
+ for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) ) { \r
+ SwapEndian_32p(&x[i]); \r
+ }\r
+ }\r
+ \r
+ // set BamAlignment 'core' and 'support' data\r
+ bAlignment.RefID = BgzfData::UnpackSignedInt(&x[0]); \r
+ bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);\r
+ \r
+ unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);\r
+ bAlignment.Bin = tempValue >> 16;\r
+ bAlignment.MapQuality = tempValue >> 8 & 0xff;\r
+ bAlignment.SupportData.QueryNameLength = tempValue & 0xff;\r
+\r
+ tempValue = BgzfData::UnpackUnsignedInt(&x[12]);\r
+ bAlignment.AlignmentFlag = tempValue >> 16;\r
+ bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;\r
+\r
+ bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);\r
+ bAlignment.MateRefID = BgzfData::UnpackSignedInt(&x[20]);\r
+ bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);\r
+ bAlignment.InsertSize = BgzfData::UnpackSignedInt(&x[28]);\r
+ \r
+ // set BamAlignment length\r
+ bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;\r
+ \r
+ // read in character data - make sure proper data size was read\r
+ bool readCharDataOK = false;\r
+ const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;\r
+ char* allCharData = (char*)calloc(sizeof(char), dataLength);\r
+ \r
+ if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) { \r
+ \r
+ // store 'allCharData' in supportData structure\r
+ bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);\r
+ \r
+ // set success flag\r
+ readCharDataOK = true;\r
+ }\r
+\r
+ free(allCharData);\r
+ return readCharDataOK;\r
+}\r
+\r
+// loads reference data from BAM file\r
+void BamReader::BamReaderPrivate::LoadReferenceData(void) {\r
+\r
+ // get number of reference sequences\r
+ char buffer[4];\r
+ mBGZF.Read(buffer, 4);\r
+ unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);\r
+ if ( IsBigEndian ) { SwapEndian_32(numberRefSeqs); }\r
+ if (numberRefSeqs == 0) { return; }\r
+ References.reserve((int)numberRefSeqs);\r
+\r
+ // iterate over all references in header\r
+ for (unsigned int i = 0; i != numberRefSeqs; ++i) {\r
+\r
+ // get length of reference name\r
+ mBGZF.Read(buffer, 4);\r
+ unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);\r
+ if ( IsBigEndian ) { SwapEndian_32(refNameLength); }\r
+ char* refName = (char*)calloc(refNameLength, 1);\r
+\r
+ // get reference name and reference sequence length\r
+ mBGZF.Read(refName, refNameLength);\r
+ mBGZF.Read(buffer, 4);\r
+ int refLength = BgzfData::UnpackSignedInt(buffer);\r
+ if ( IsBigEndian ) { SwapEndian_32(refLength); }\r
+\r
+ // store data for reference\r
+ RefData aReference;\r
+ aReference.RefName = (string)((const char*)refName);\r
+ aReference.RefLength = refLength;\r
+ References.push_back(aReference);\r
+\r
+ // clean up calloc-ed temp variable\r
+ free(refName);\r
+ }\r
+}\r
+\r
+// opens BAM file (and index)\r
+bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename) {\r
+\r
+ Filename = filename;\r
+ IndexFilename = indexFilename;\r
+\r
+ // open the BGZF file for reading, return false on failure\r
+ if ( !mBGZF.Open(filename, "rb") ) \r
+ return false;\r
+ \r
+ // retrieve header text & reference data\r
+ LoadHeaderData();\r
+ LoadReferenceData();\r
+\r
+ // store file offset of first alignment\r
+ AlignmentsBeginOffset = mBGZF.Tell();\r
+\r
+ // open index file & load index data (if exists)\r
+ if ( !IndexFilename.empty() )\r
+ LoadIndex();\r
+ \r
+ // return success\r
+ return true;\r
+}\r
+\r
+// returns BAM file pointer to beginning of alignment data\r
+bool BamReader::BamReaderPrivate::Rewind(void) {\r
+ \r
+ // rewind to first alignment\r
+ if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;\r
+ \r
+ // retrieve first alignment data\r
+ BamAlignment al;\r
+ if ( !LoadNextAlignment(al) ) return false;\r
+ \r
+ // reset default region info using first alignment in file\r
+ Region.LeftRefID = al.RefID;\r
+ Region.LeftPosition = al.Position;\r
+ Region.RightRefID = -1;\r
+ Region.RightPosition = -1;\r
+ IsLeftBoundSpecified = false;\r
+ IsRightBoundSpecified = false; \r
+\r
+ // rewind back to before first alignment\r
+ // return success/fail of seek\r
+ return mBGZF.Seek(AlignmentsBeginOffset);\r
+}\r
+\r
+// sets a region of interest (with left & right bound reference/position)\r
+// attempts a Jump() to left bound as well\r
+// returns success/failure of Jump()\r
+bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {\r
+ \r
+ // save region of interest\r
+ Region = region;\r
+ \r
+ // set flags\r
+ if ( region.LeftRefID >= 0 && region.LeftPosition >= 0 ) \r
+ IsLeftBoundSpecified = true;\r
+ if ( region.RightRefID >= 0 && region.RightPosition >= 0 ) \r
+ IsRightBoundSpecified = true;\r
+ \r
+ // attempt jump to beginning of region, return success/fail of Jump()\r
+ return Jump( Region.LeftRefID, Region.LeftPosition );\r
+}\r
projects / bamtools.git / blobdiff