-// ***************************************************************************\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
- fprintf(stderr, "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
- fprintf(stderr, "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
- fprintf(stderr, "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
- fprintf(stderr, "Could not read header type\n");\r
- exit(1);\r
- }\r
-\r
- if (strncmp(buffer, "BAM\001", 4)) {\r
- fprintf(stderr, "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
- fprintf(stderr, "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
+// ***************************************************************************
+// BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg
+// Marth Lab, Department of Biology, Boston College
+// ---------------------------------------------------------------------------
+// Last modified: 10 October 2011 (DB)
+// ---------------------------------------------------------------------------
+// Provides read access to BAM files.
+// ***************************************************************************
+
+#include <api/BamReader.h>
+#include <api/internal/BamReader_p.h>
+using namespace BamTools;
+using namespace BamTools::Internal;
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <string>
+#include <vector>
+using namespace std;
+
+/*! \class BamTools::BamReader
+ \brief Provides read access to BAM files.
+*/
+
+/*! \fn BamReader::BamReader(void)
+ \brief constructor
+*/
+BamReader::BamReader(void)
+ : d(new BamReaderPrivate(this))
+{ }
+
+/*! \fn BamReader::~BamReader(void)
+ \brief destructor
+*/
+BamReader::~BamReader(void) {
+ delete d;
+ d = 0;
+}
+
+/*! \fn bool BamReader::Close(void)
+ \brief Closes the current BAM file.
+
+ Also clears out all header and reference data.
+
+ \return \c true if file closed OK
+ \sa IsOpen(), Open()
+*/
+bool BamReader::Close(void) {
+ return d->Close();
+}
+
+/*! \fn bool BamReader::CreateIndex(const BamIndex::IndexType& type)
+ \brief Creates an index file for current BAM file.
+
+ \param[in] type file format to create, see BamIndex::IndexType for available formats
+ \return \c true if index created OK
+ \sa LocateIndex(), OpenIndex()
+*/
+bool BamReader::CreateIndex(const BamIndex::IndexType& type) {
+ return d->CreateIndex(type);
+}
+
+/*! \fn std::string BamReader::GetErrorString(void) const
+ \brief Returns a human-readable description of the last error that occurred
+
+ This method allows elimination of STDERR pollution. Developers of client code
+ may choose how the messages are displayed to the user, if at all.
+
+ \return error description
+*/
+string BamReader::GetErrorString(void) const {
+ return d->GetErrorString();
+}
+
+/*! \fn const std::string BamReader::GetFilename(void) const
+ \brief Returns name of current BAM file.
+
+ Retrieved filename will contain whatever was passed via Open().
+ If you need full directory paths here, be sure to include them
+ when you open the BAM file.
+
+ \returns name of open BAM file. If no file is open, returns an empty string.
+ \sa IsOpen()
+*/
+const std::string BamReader::GetFilename(void) const {
+ return d->Filename();
+}
+
+/*! \fn SamHeader BamReader::GetHeader(void) const
+ \brief Returns SAM header data.
+
+ Header data is wrapped in a SamHeader object that can be conveniently queried & modified.
+
+ \note Modifying the retrieved SamHeader object does NOT affect the
+ current BAM file. This file has been opened in a read-only mode.
+ However, your modified SamHeader object can be used in conjunction with
+ BamWriter to generate a new BAM file with the appropriate header information.
+
+ \returns header data object
+ \sa GetHeaderText()
+*/
+SamHeader BamReader::GetHeader(void) const {
+ return d->GetSamHeader();
+}
+
+/*! \fn std::string BamReader::GetHeaderText(void) const
+ \brief Returns SAM header data, as SAM-formatted text.
+
+ \note Modifying the retrieved text does NOT affect the current
+ BAM file. This file has been opened in a read-only mode. However,
+ your modified header text can be used in conjunction with BamWriter
+ to generate a new BAM file with the appropriate header information.
+
+ \returns SAM-formatted header text
+ \sa GetHeader()
+*/
+std::string BamReader::GetHeaderText(void) const {
+ return d->GetHeaderText();
+}
+
+/*! \fn bool BamReader::GetNextAlignment(BamAlignment& alignment)
+ \brief Retrieves next available alignment.
+
+ Attempts to read the next alignment record from BAM file, and checks to see
+ if it overlaps the current region. If no region is currently set, then the
+ next alignment available is always considered valid.
+
+ If a region has been set, via Jump() or SetRegion(), an alignment is only
+ considered valid if it overlaps the region. If the actual 'next' alignment record
+ in the BAM file does not overlap this region, then this function will read sequentially
+ through the file until the next alignment that overlaps this region is found.
+ Once the region has been exhausted (i.e. the next alignment loaded is beyond the region),
+ the function aborts and returns \c false. In this case, there is no point to continue
+ reading, assuming properly sorted alignments.
+
+ This function fully populates all of the alignment's available data fields,
+ including the string data fields (read name, bases, qualities, tags, filename).
+ If only positional data (refID, position, CIGAR ops, alignment flags, etc.)
+ are required, consider using GetNextAlignmentCore() for a significant
+ performance boost.
+
+ \param[out] alignment destination for alignment record data
+ \returns \c true if a valid alignment was found
+*/
+bool BamReader::GetNextAlignment(BamAlignment& alignment) {
+ return d->GetNextAlignment(alignment);
+}
+
+/*! \fn bool BamReader::GetNextAlignmentCore(BamAlignment& alignment)
+ \brief Retrieves next available alignment, without populating the alignment's string data fields.
+
+ Equivalent to GetNextAlignment() with respect to what is a valid overlapping alignment.
+
+ However, this method does NOT populate the alignment's string data fields
+ (read name, bases, qualities, tags, filename). This provides a boost in speed
+ when these fields are not required for every alignment. These fields can be
+ populated 'lazily' (as needed) by calling BamAlignment::BuildCharData() later.
+
+ \param[out] alignment destination for alignment record data
+ \returns \c true if a valid alignment was found
+ \sa SetRegion()
+*/
+bool BamReader::GetNextAlignmentCore(BamAlignment& alignment) {
+ return d->GetNextAlignmentCore(alignment);
+}
+
+/*! \fn int BamReader::GetReferenceCount(void) const
+ \brief Returns number of reference sequences.
+*/
+int BamReader::GetReferenceCount(void) const {
+ return d->GetReferenceCount();
+}
+
+/*! \fn const RefVector& BamReader::GetReferenceData(void) const
+ \brief Returns all reference sequence entries.
+ \sa RefData
+*/
+const RefVector& BamReader::GetReferenceData(void) const {
+ return d->GetReferenceData();
+}
+
+/*! \fn int BamReader::GetReferenceID(const std::string& refName) const
+ \brief Returns the ID of the reference with this name.
+
+ If \a refName is not found, returns -1.
+
+ \param[in] refName name of reference to look up
+*/
+int BamReader::GetReferenceID(const std::string& refName) const {
+ return d->GetReferenceID(refName);
+}
+
+/*! \fn bool BamReader::HasIndex(void) const
+ \brief Returns \c true if index data is available.
+*/
+bool BamReader::HasIndex(void) const {
+ return d->HasIndex();
+}
+
+/*! \fn bool BamReader::IsOpen(void) const
+ \brief Returns \c true if a BAM file is open for reading.
+*/
+bool BamReader::IsOpen(void) const {
+ return d->IsOpen();
+}
+
+/*! \fn bool BamReader::Jump(int refID, int position)
+ \brief Performs a random-access jump within BAM file.
+
+ This is a convenience method, equivalent to calling SetRegion()
+ with only a left boundary specified.
+
+ \param[in] refID left-bound reference ID
+ \param[in] position left-bound position
+
+ \returns \c true if jump was successful
+ \sa HasIndex()
+*/
+bool BamReader::Jump(int refID, int position) {
+ return d->SetRegion( BamRegion(refID, position) );
+}
+
+/*! \fn bool BamReader::LocateIndex(const BamIndex::IndexType& preferredType)
+ \brief Looks in BAM file's directory for a matching index file.
+
+ Use this function when you need an index file, and perhaps have a
+ preferred index format, but do not depend heavily on which format
+ actually gets loaded at runtime.
+
+ This function will defer to your \a preferredType whenever possible.
+ However, if an index file of \a preferredType can not be found, then
+ it will look for any other index file that corresponds to this BAM file.
+
+ If you want precise control over which index file is loaded, use OpenIndex()
+ with the desired index filename. If that function returns false, you can use
+ CreateIndex() to then build an index of the exact requested format.
+
+ \param[in] preferredType desired index file format, see BamIndex::IndexType for available formats
+
+ \returns \c true if (any) index file could be found
+*/
+bool BamReader::LocateIndex(const BamIndex::IndexType& preferredType) {
+ return d->LocateIndex(preferredType);
+}
+
+/*! \fn bool BamReader::Open(const std::string& filename)
+ \brief Opens a BAM file.
+
+ If BamReader is already opened on another file, this function closes
+ that file, then attempts to open requested \a filename.
+
+ \param[in] filename name of BAM file to open
+
+ \returns \c true if BAM file was opened successfully
+ \sa Close(), IsOpen(), OpenIndex()
+*/
+bool BamReader::Open(const std::string& filename) {
+ return d->Open(filename);
+}
+
+/*! \fn bool BamReader::OpenIndex(const std::string& indexFilename)
+ \brief Opens a BAM index file.
+
+ \param[in] indexFilename name of BAM index file to open
+
+ \returns \c true if BAM index file was opened & data loaded successfully
+ \sa LocateIndex(), Open(), SetIndex()
+*/
+bool BamReader::OpenIndex(const std::string& indexFilename) {
+ return d->OpenIndex(indexFilename);
+}
+
+/*! \fn bool BamReader::Rewind(void)
+ \brief Returns the internal file pointer to the first alignment record.
+
+ Useful for performing multiple sequential passes through a BAM file.
+ Calling this function clears any prior region that may have been set.
+
+ \note This function sets the file pointer to first alignment record
+ in the BAM file, NOT the beginning of the file.
+
+ \returns \c true if rewind operation was successful
+ \sa Jump(), SetRegion()
+*/
+bool BamReader::Rewind(void) {
+ return d->Rewind();
+}
+
+/*! \fn void BamReader::SetIndex(BamIndex* index)
+ \brief Sets a custom BamIndex on this reader.
+
+ Only necessary for custom BamIndex subclasses. Most clients should
+ never have to use this function.
+
+ Example:
+ \code
+ BamReader reader;
+ reader.SetIndex(new MyCustomBamIndex);
+ \endcode
+
+ \note BamReader takes ownership of \a index - i.e. the BamReader will
+ take care of deleting it when the reader is destructed, when the current
+ BAM file is closed, or when a new index is requested.
+
+ \param[in] index custom BamIndex subclass created by client
+ \sa CreateIndex(), LocateIndex(), OpenIndex()
+*/
+void BamReader::SetIndex(BamIndex* index) {
+ d->SetIndex(index);
+}
+
+/*! \fn void BamReader::SetIndexCacheMode(const BamIndex::IndexCacheMode& mode)
+ \brief Changes the caching behavior of the index data.
+
+ Default mode is BamIndex::LimitedIndexCaching.
+
+ \param[in] mode desired cache mode for index, see BamIndex::IndexCacheMode for
+ description of the available cache modes
+ \sa HasIndex()
+*/
+void BamReader::SetIndexCacheMode(const BamIndex::IndexCacheMode& mode) {
+ d->SetIndexCacheMode(mode);
+}
+
+/*! \fn bool BamReader::SetRegion(const BamRegion& region)
+ \brief Sets a target region of interest
+
+ Requires that index data be available. Attempts a random-access
+ jump in the BAM file, near \a region left boundary position.
+
+ Subsequent calls to GetNextAlignment() or GetNextAlignmentCore()
+ will only return \c true when alignments can be found that overlap
+ this \a region.
+
+ A \a region with no right boundary is considered open-ended, meaning
+ that all alignments that lie downstream of the left boundary are
+ considered valid, continuing to the end of the BAM file.
+
+ \warning BamRegion now represents a zero-based, HALF-OPEN interval.
+ In previous versions of BamTools (0.x & 1.x) all intervals were treated
+ as zero-based, CLOSED.
+
+ \param[in] region desired region-of-interest to activate
+
+ \returns \c true if reader was able to jump successfully to the region's left boundary
+ \sa HasIndex(), Jump()
+*/
+bool BamReader::SetRegion(const BamRegion& region) {
+ return d->SetRegion(region);
+}
+
+/*! \fn bool BamReader::SetRegion(const int& leftRefID,
+ const int& leftPosition,
+ const int& rightRefID,
+ const int& rightPosition)
+ \brief Sets a target region of interest.
+
+ This is an overloaded function.
+
+ \warning This function expects a zero-based, HALF-OPEN interval.
+ In previous versions of BamTools (0.x & 1.x) all intervals were treated
+ as zero-based, CLOSED.
+
+ \param[in] leftRefID referenceID of region's left boundary
+ \param[in] leftPosition position of region's left boundary
+ \param[in] rightRefID reference ID of region's right boundary
+ \param[in] rightPosition position of region's right boundary
+
+ \returns \c true if reader was able to jump successfully to the region's left boundary
+ \sa HasIndex(), Jump()
+*/
+bool BamReader::SetRegion(const int& leftRefID,
+ const int& leftBound,
+ const int& rightRefID,
+ const int& rightBound)
+{
+ return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );
+}
projects / bamtools.git / blobdiff