X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fapi%2FBamMultiReader.cpp;h=57c826daa80ccd5368bc0cde3c7254d1f417f496;hb=574a2bfb36f7107529e7ccda0f75e70a493460e5;hp=6de97dc9e7f15411e887ea1c97a62d41fe1ed9d5;hpb=164366c9fd23717e4b7279518eecaa32baed502c;p=bamtools.git diff --git a/src/api/BamMultiReader.cpp b/src/api/BamMultiReader.cpp index 6de97dc..57c826d 100644 --- a/src/api/BamMultiReader.cpp +++ b/src/api/BamMultiReader.cpp @@ -1,467 +1,417 @@ // *************************************************************************** -// BamMultiReader.cpp (c) 2010 Erik Garrison +// BamMultiReader.cpp (c) 2010 Erik Garrison, Derek Barnett // Marth Lab, Department of Biology, Boston College -// All rights reserved. // --------------------------------------------------------------------------- -// Last modified: 3 September 2010 (DB) +// Last modified: 14 January 2013 (DB) // --------------------------------------------------------------------------- -// Uses BGZF routines were adapted from the bgzf.c code developed at the Broad -// Institute. -// --------------------------------------------------------------------------- -// Functionality for simultaneously reading multiple BAM files. +// Convenience class for reading multiple BAM files. // // This functionality allows applications to work on very large sets of files // without requiring intermediate merge, sort, and index steps for each file -// subset. It also improves the performance of our merge system as it +// subset. It also improves the performance of our merge system as it // precludes the need to sort merged files. // *************************************************************************** -// C++ includes -#include -#include -#include -#include -#include +#include "api/BamMultiReader.h" +#include "api/internal/bam/BamMultiReader_p.h" +using namespace BamTools; + #include #include - -// BamTools includes -#include "BGZF.h" -#include "BamMultiReader.h" -using namespace BamTools; using namespace std; -namespace BamTools { +/*! \class BamTools::BamMultiReader + \brief Convenience class for reading multiple BAM files. +*/ + +/*! \enum BamMultiReader::MergeOrder + \brief A description of the enum type. +*/ +/*! \var BamMultiReader::MergeOrder BamMultiReader::MergeByCoordinate + \brief The description of the first enum value. +*/ +/*! \var BamMultiReader::MergeOrder BamMultiReader::MergeByName + \brief BAM files are +*/ -bool FileExists(const string& filename) { - ifstream f(filename.c_str(), ifstream::in); - return !f.fail(); -} - -} // namespace BamTools -// ----------------------------------------------------- -// BamMultiReader implementation -// ----------------------------------------------------- -// constructor +/*! \fn BamMultiReader::BamMultiReader(void) + \brief constructor +*/ BamMultiReader::BamMultiReader(void) - : CurrentRefID(0) - , CurrentLeft(0) + : d(new Internal::BamMultiReaderPrivate) { } -// destructor +/*! \fn BamMultiReader::~BamMultiReader(void) + \brief destructor +*/ BamMultiReader::~BamMultiReader(void) { - Close(); + delete d; + d = 0; } -// close the BAM files -void BamMultiReader::Close(void) { - - // close all BAM readers and clean up pointers - vector >::iterator readerIter = readers.begin(); - vector >::iterator readerEnd = readers.end(); - for ( ; readerIter != readerEnd; ++readerIter) { - - BamReader* reader = (*readerIter).first; - BamAlignment* alignment = (*readerIter).second; - - // close the reader - if ( reader) reader->Close(); - - // delete reader pointer - delete reader; - reader = 0; - - // delete alignment pointer - delete alignment; - alignment = 0; - } - - // clear out the container - readers.clear(); -} +/*! \fn void BamMultiReader::Close(void) + \brief Closes all open BAM files. + + Also clears out all header and reference data. -// saves index data to BAM index files (".bai"/".bti") where necessary, returns success/fail -bool BamMultiReader::CreateIndexes(bool useDefaultIndex) { - bool result = true; - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* reader = it->first; - result &= reader->CreateIndex(useDefaultIndex); - } - return result; + \sa CloseFile(), IsOpen(), Open(), BamReader::Close() +*/ +bool BamMultiReader::Close(void) { + return d->Close(); } -// for debugging -void BamMultiReader::DumpAlignmentIndex(void) { - for (AlignmentIndex::const_iterator it = alignments.begin(); it != alignments.end(); ++it) { - cerr << it->first.first << ":" << it->first.second << " " << it->second.first->GetFilename() << endl; - } +/*! \fn void BamMultiReader::CloseFile(const std::string& filename) + \brief Closes requested BAM file. + + Leaves any other file(s) open, along with header and reference data. + + \param[in] filename name of specific BAM file to close + + \sa Close(), IsOpen(), Open(), BamReader::Close() +*/ +bool BamMultiReader::CloseFile(const std::string& filename) { + return d->CloseFile(filename); } -// makes a virtual, unified header for all the bam files in the multireader -const string BamMultiReader::GetHeaderText(void) const { - - string mergedHeader = ""; - map readGroups; - - // foreach extraction entry (each BAM file) - for (vector >::const_iterator rs = readers.begin(); rs != readers.end(); ++rs) { - - map currentFileReadGroups; - - BamReader* reader = rs->first; - - stringstream header(reader->GetHeaderText()); - vector lines; - string item; - while (getline(header, item)) - lines.push_back(item); - - for (vector::const_iterator it = lines.begin(); it != lines.end(); ++it) { - - // get next line from header, skip if empty - string headerLine = *it; - if ( headerLine.empty() ) { continue; } - - // if first file, save HD & SQ entries - if ( rs == readers.begin() ) { - if ( headerLine.find("@HD") == 0 || headerLine.find("@SQ") == 0) { - mergedHeader.append(headerLine.c_str()); - mergedHeader.append(1, '\n'); - } - } - - // (for all files) append RG entries if they are unique - if ( headerLine.find("@RG") == 0 ) { - stringstream headerLineSs(headerLine); - string part, readGroupPart, readGroup; - while(std::getline(headerLineSs, part, '\t')) { - stringstream partSs(part); - string subtag; - std::getline(partSs, subtag, ':'); - if (subtag == "ID") { - std::getline(partSs, readGroup, ':'); - break; - } - } - if (readGroups.find(readGroup) == readGroups.end()) { // prevents duplicate @RG entries - mergedHeader.append(headerLine.c_str() ); - mergedHeader.append(1, '\n'); - readGroups[readGroup] = true; - currentFileReadGroups[readGroup] = true; - } else { - // warn iff we are reading one file and discover duplicated @RG tags in the header - // otherwise, we emit no warning, as we might be merging multiple BAM files with identical @RG tags - if (currentFileReadGroups.find(readGroup) != currentFileReadGroups.end()) { - cerr << "WARNING: duplicate @RG tag " << readGroup - << " entry in header of " << reader->GetFilename() << endl; - } - } - } - } - } - - // return merged header text - return mergedHeader; +/*! \fn bool BamMultiReader::CreateIndexes(const BamIndex::IndexType& type) + \brief Creates index files for the current BAM files. + + \param[in] type file format to create, see BamIndex::IndexType for available formats + \return \c true if index files created OK + \sa LocateIndexes(), OpenIndexes(), BamReader::CreateIndex() +*/ +bool BamMultiReader::CreateIndexes(const BamIndex::IndexType& type) { + return d->CreateIndexes(type); } -// get next alignment among all files -bool BamMultiReader::GetNextAlignment(BamAlignment& nextAlignment) { +/*! \fn const std::vector BamMultiReader::Filenames(void) const + \brief Returns list of filenames for all open BAM files. - // bail out if we are at EOF in all files, means no more alignments to process - if (!HasOpenReaders()) - return false; + Retrieved filenames will contain whatever was passed via Open(). + If you need full directory paths here, be sure to include them + when you open the BAM files. - // when all alignments have stepped into a new target sequence, update our - // current reference sequence id - UpdateReferenceID(); + \returns names of open BAM files. If no files are open, returns an empty vector. + \sa IsOpen(), BamReader::GetFilename() +*/ +const std::vector BamMultiReader::Filenames(void) const { + return d->Filenames(); +} - // our lowest alignment and reader will be at the front of our alignment index - BamAlignment* alignment = alignments.begin()->second.second; - BamReader* reader = alignments.begin()->second.first; +/*! \fn std::string BamMultiReader::GetErrorString(void) const + \brief Returns a human-readable description of the last error that occurred - // now that we have the lowest alignment in the set, save it by copy to our argument - nextAlignment = BamAlignment(*alignment); + This method allows elimination of STDERR pollution. Developers of client code + may choose how the messages are displayed to the user, if at all. - // remove this alignment index entry from our alignment index - alignments.erase(alignments.begin()); + \return error description +*/ +std::string BamMultiReader::GetErrorString(void) const { + return d->GetErrorString(); +} - // and add another entry if we can get another alignment from the reader - if (reader->GetNextAlignment(*alignment)) { - alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position), - make_pair(reader, alignment))); - } else { // do nothing - //cerr << "reached end of file " << lowestReader->GetFilename() << endl; - } +/*! \fn SamHeader BamMultiReader::GetHeader(void) const + \brief Returns unified SAM-format header for all files - return true; + \note Modifying the retrieved text does NOT affect the current + BAM files. These files have 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 header data wrapped in SamHeader object + \sa GetHeaderText(), BamReader::GetHeader() +*/ +SamHeader BamMultiReader::GetHeader(void) const { + return d->GetHeader(); } -// get next alignment among all files without parsing character data from alignments -bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) { +/*! \fn std::string BamMultiReader::GetHeaderText(void) const + \brief Returns unified SAM-format header text for all files - // bail out if we are at EOF in all files, means no more alignments to process - if (!HasOpenReaders()) - return false; + \note Modifying the retrieved text does NOT affect the current + BAM files. These files have 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. - // when all alignments have stepped into a new target sequence, update our - // current reference sequence id - UpdateReferenceID(); + \returns SAM-formatted header text + \sa GetHeader(), BamReader::GetHeaderText() +*/ +std::string BamMultiReader::GetHeaderText(void) const { + return d->GetHeaderText(); +} - // our lowest alignment and reader will be at the front of our alignment index - BamAlignment* alignment = alignments.begin()->second.second; - BamReader* reader = alignments.begin()->second.first; +/*! \fn BamMultiReader::MergeOrder BamMultiReader::GetMergeOrder(void) const + \brief Returns curent merge order strategy. - // now that we have the lowest alignment in the set, save it by copy to our argument - nextAlignment = BamAlignment(*alignment); - //memcpy(&nextAlignment, alignment, sizeof(BamAlignment)); + \returns current merge order enum value + \sa BamMultiReader::MergeOrder, SetExplicitMergeOrder() +*/ +BamMultiReader::MergeOrder BamMultiReader::GetMergeOrder(void) const { + return d->GetMergeOrder(); +} - // remove this alignment index entry from our alignment index - alignments.erase(alignments.begin()); +/*! \fn bool BamMultiReader::GetNextAlignment(BamAlignment& alignment) + \brief Retrieves next available alignment. - // and add another entry if we can get another alignment from the reader - if (reader->GetNextAlignmentCore(*alignment)) { - alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position), - make_pair(reader, alignment))); - } else { // do nothing - //cerr << "reached end of file " << lowestReader->GetFilename() << endl; - } + Equivalent to BamReader::GetNextAlignment() with respect to what is a valid + overlapping alignment and what data gets populated. - return true; + This method takes care of determining which alignment actually is 'next' + across multiple files, depending on their sort order. + \param[out] alignment destination for alignment record data + \returns \c true if a valid alignment was found + \sa GetNextAlignmentCore(), SetExplicitMergeOrder(), SetRegion(), BamReader::GetNextAlignment() +*/ +bool BamMultiReader::GetNextAlignment(BamAlignment& nextAlignment) { + return d->GetNextAlignment(nextAlignment); } -// --------------------------------------------------------------------------------------- -// -// NB: The following GetReferenceX() functions assume that we have identical -// references for all BAM files. We enforce this by invoking the above -// validation function (ValidateReaders) to verify that our reference data -// is the same across all files on Open, so we will not encounter a situation -// in which there is a mismatch and we are still live. -// -// --------------------------------------------------------------------------------------- +/*! \fn bool BamMultiReader::GetNextAlignmentCore(BamAlignment& alignment) + \brief Retrieves next available alignment. + + Equivalent to BamReader::GetNextAlignmentCore() with respect to what is a valid + overlapping alignment and what data gets populated. -// returns the number of reference sequences -const int BamMultiReader::GetReferenceCount(void) const { - return readers.front().first->GetReferenceCount(); + This method takes care of determining which alignment actually is 'next' + across multiple files, depending on their sort order. + + \param[out] alignment destination for alignment record data + \returns \c true if a valid alignment was found + \sa GetNextAlignment(), SetExplicitMergeOrder(), SetRegion(), BamReader::GetNextAlignmentCore() +*/ +bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) { + return d->GetNextAlignmentCore(nextAlignment); } -// returns vector of reference objects -const BamTools::RefVector BamMultiReader::GetReferenceData(void) const { - return readers.front().first->GetReferenceData(); +/*! \fn int BamMultiReader::GetReferenceCount(void) const + \brief Returns number of reference sequences. + \sa BamReader::GetReferenceCount() +*/ +int BamMultiReader::GetReferenceCount(void) const { + return d->GetReferenceCount(); } -// returns refID from reference name -const int BamMultiReader::GetReferenceID(const string& refName) const { - return readers.front().first->GetReferenceID(refName); +/*! \fn const RefVector& BamMultiReader::GetReferenceData(void) const + \brief Returns all reference sequence entries. + \sa RefData, BamReader::GetReferenceData() +*/ +const BamTools::RefVector BamMultiReader::GetReferenceData(void) const { + return d->GetReferenceData(); } -// --------------------------------------------------------------------------------------- +/*! \fn int BamMultiReader::GetReferenceID(const std::string& refName) const + \brief Returns the ID of the reference with this name. -// checks if any readers still have alignments -bool BamMultiReader::HasOpenReaders() { - return alignments.size() > 0; -} + If \a refName is not found, returns -1. -// returns whether underlying BAM readers ALL have an index loaded -// this is useful to indicate whether Jump() or SetRegion() are possible -bool BamMultiReader::IsIndexLoaded(void) const { - bool ok = true; - vector >::const_iterator readerIter = readers.begin(); - vector >::const_iterator readerEnd = readers.end(); - for ( ; readerIter != readerEnd; ++readerIter ) { - const BamReader* reader = (*readerIter).first; - if ( reader ) ok &= reader->IsIndexLoaded(); - } - return ok; + \param[in] refName name of reference to look up + \sa BamReader::GetReferenceID() +*/ +int BamMultiReader::GetReferenceID(const std::string& refName) const { + return d->GetReferenceID(refName); } -// jumps to specified region(refID, leftBound) in BAM files, returns success/fail -bool BamMultiReader::Jump(int refID, int position) { +/*! \fn bool BamMultiReader::HasIndexes(void) const + \brief Returns \c true if all BAM files have index data available. + \sa BamReader::HasIndex() +*/ +bool BamMultiReader::HasIndexes(void) const { + return d->HasIndexes(); +} - //if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) { - CurrentRefID = refID; - CurrentLeft = position; - - bool result = true; - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* reader = it->first; - result &= reader->Jump(refID, position); - if (!result) { - cerr << "ERROR: could not jump " << reader->GetFilename() << " to " << refID << ":" << position << endl; - exit(1); - } - } - if (result) UpdateAlignments(); - return result; +/*! \fn bool BamMultiReader::HasOpenReaders(void) const + \brief Returns \c true if there are any open BAM files. +*/ +bool BamMultiReader::HasOpenReaders(void) const { + return d->HasOpenReaders(); } -// opens BAM files -bool BamMultiReader::Open(const vector filenames, bool openIndexes, bool coreMode, bool useDefaultIndex) { - - // for filename in filenames - fileNames = filenames; // save filenames in our multireader - for (vector::const_iterator it = filenames.begin(); it != filenames.end(); ++it) { - - const string filename = *it; - BamReader* reader = new BamReader; - - bool openedOK = true; - if (openIndexes) { - - const string defaultIndexFilename = filename + ".bai"; - const string bamToolsIndexFilename = filename + ".bti"; - - // if default BAM index requested and one exists - if ( useDefaultIndex && FileExists(defaultIndexFilename) ) - openedOK = reader->Open(filename, defaultIndexFilename); - - // else see if BamTools index exists - else if ( FileExists(bamToolsIndexFilename) ) - openedOK = reader->Open(filename, bamToolsIndexFilename); - - // else none exist... just open without - else - openedOK = reader->Open(filename); - } - - // ignoring index file(s) - else openedOK = reader->Open(filename); - - // if file opened ok, check that it can be read - if ( openedOK ) { - - bool fileOK = true; - BamAlignment* alignment = new BamAlignment; - fileOK &= ( coreMode ? reader->GetNextAlignmentCore(*alignment) : reader->GetNextAlignment(*alignment) ); - - if (fileOK) { - readers.push_back(make_pair(reader, alignment)); // store pointers to our readers for cleanup - alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position), - make_pair(reader, alignment))); - } else { - cerr << "WARNING: could not read first alignment in " << filename << ", ignoring file" << endl; - // if only file available & could not be read, return failure - if ( filenames.size() == 1 ) return false; - } - } - - // TODO; any more error handling when openedOKvis false ?? - else - return false; - } - - // files opened ok, at least one alignment could be read, - // now need to check that all files use same reference data - ValidateReaders(); - return true; +/*! \fn bool BamMultiReader::Jump(int refID, int position) + \brief Performs a random-access jump within current BAM files. + + This is a convenience method, equivalent to calling SetRegion() + with only a left boundary specified. + + \param[in] refID ID of reference to jump to + \param[in] position (0-based) left boundary + + \returns \c true if jump was successful + \sa HasIndex(), BamReader::Jump() +*/ + +bool BamMultiReader::Jump(int refID, int position) { + return d->Jump(refID, position); } -void BamMultiReader::PrintFilenames(void) { - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* reader = it->first; - cout << reader->GetFilename() << endl; - } +/*! \fn bool BamMultiReader::LocateIndexes(const BamIndex::IndexType& preferredType) + \brief Looks for index files that match current BAM files. + + Use this function when you need index files, and perhaps have a + preferred index format, but do not depend heavily on which indexes + actually get loaded at runtime. + + For each BAM file, 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 matches + that BAM file. + + An example case would look this: + \code + BamMultiReader reader; + + // do setup... + + // ensure that all files have an index + if ( !reader.LocateIndexes() ) // opens any existing index files that match our BAM files + reader.CreateIndexes(); // creates index files for any BAM files that still lack one + + // do interesting stuff using random-access... + + \endcode + + If you want precise control over which index files are loaded, use OpenIndexes() + with the desired index filenames. If that function returns false, you can use + CreateIndexes() to then build index files of the exact requested format. + + \param[in] preferredType desired index file format, see BamIndex::IndexType for available formats + \returns \c true if index files could be found for \b ALL open BAM files + \sa BamReader::LocateIndex() +*/ +bool BamMultiReader::LocateIndexes(const BamIndex::IndexType& preferredType) { + return d->LocateIndexes(preferredType); } -// returns BAM file pointers to beginning of alignment data -bool BamMultiReader::Rewind(void) { - bool result = true; - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* reader = it->first; - result &= reader->Rewind(); - } - return result; +/*! \fn bool BamMultiReader::Open(const std::vector& filenames) + \brief Opens BAM files. + + \note Opening BAM files will invalidate any current region set on the multireader. + All file pointers will be returned to the beginning of the alignment data. Follow + this with Jump() or SetRegion() to establish a region of interest. + + \param[in] filenames list of BAM filenames to open + \returns \c true if BAM files were opened successfully + \sa Close(), HasOpenReaders(), OpenFile(), OpenIndexes(), BamReader::Open() +*/ +bool BamMultiReader::Open(const std::vector& filenames) { + return d->Open(filenames); } -bool BamMultiReader::SetRegion(const int& leftRefID, const int& leftPosition, const int& rightRefID, const int& rightPosition) { - BamRegion region(leftRefID, leftPosition, rightRefID, rightPosition); - return SetRegion(region); +/*! \fn bool BamMultiReader::OpenFile(const std::string& filename) + \brief Opens a single BAM file. + + Adds another BAM file to multireader "on-the-fly". + + \note Opening a BAM file will invalidate any current region set on the multireader. + All file pointers will be returned to the beginning of the alignment data. Follow + this with Jump() or SetRegion() to establish a region of interest. + + \param[in] filename BAM filename to open + \returns \c true if BAM file was opened successfully + \sa Close(), HasOpenReaders(), Open(), OpenIndexes(), BamReader::Open() +*/ +bool BamMultiReader::OpenFile(const std::string& filename) { + return d->OpenFile(filename); } -bool BamMultiReader::SetRegion(const BamRegion& region) { +/*! \fn bool BamMultiReader::OpenIndexes(const std::vector& indexFilenames) + \brief Opens index files for current BAM files. + + \note Currently assumes that index filenames match the order (and number) of + BAM files passed to Open(). - Region = region; + \param[in] indexFilenames list of BAM index file names + \returns \c true if BAM index file was opened & data loaded successfully + \sa LocateIndex(), Open(), SetIndex(), BamReader::OpenIndex() +*/ +bool BamMultiReader::OpenIndexes(const std::vector& indexFilenames) { + return d->OpenIndexes(indexFilenames); +} - // NB: While it may make sense to track readers in which we can - // successfully SetRegion, In practice a failure of SetRegion means "no - // alignments here." It makes sense to simply accept the failure, - // UpdateAlignments(), and continue. +/*! \fn bool BamMultiReader::Rewind(void) + \brief Returns the internal file pointers to the beginning of alignment records. - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - it->first->SetRegion(region); - } + Useful for performing multiple sequential passes through BAM files. + Calling this function clears any prior region that may have been set. - UpdateAlignments(); - return true; + \returns \c true if rewind operation was successful + \sa Jump(), SetRegion(), BamReader::Rewind() +*/ +bool BamMultiReader::Rewind(void) { + return d->Rewind(); } -void BamMultiReader::UpdateAlignments(void) { - // Update Alignments - alignments.clear(); - for (vector >::iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* br = it->first; - BamAlignment* ba = it->second; - if (br->GetNextAlignment(*ba)) { - alignments.insert(make_pair(make_pair(ba->RefID, ba->Position), - make_pair(br, ba))); - } else { - // assume BamReader end of region / EOF - } - } +/*! \fn void BamMultiReader::SetExplicitMergeOrder(BamMultiReader::MergeOrder order) + \brief Sets an explicit merge order, regardless of the BAM files' SO header tag. + + The default behavior of the BamMultiReader is to check the SO tag in the BAM files' + SAM header text to determine the merge strategy". The merge strategy is used to + determine from which BAM file the next alignment should come when either + GetNextAlignment() or GetNextAlignmentCore() are called. If files share a + 'coordinate' or 'queryname' value for this tag, then the merge strategy is + selected accordingly. If any of them do not match, or if any fileis marked as + 'unsorted', then the merge strategy is simply a round-robin. + + This method allows client code to explicitly override the lookup behavior. This + method can be useful when you know, for example, that your BAM files are sorted + by coordinate but upstream processes did not set the header tag properly. + + \note This method should \bold not be called while reading alignments via + GetNextAlignment() or GetNextAlignmentCore(). For proper results, you should + call this method before (or immediately after) opening files, rewinding, + jumping, etc. but \bold not once alignment fetching has started. There is + nothing in the API to prevent you from doing so, but the results may be + unexpected. + + \sa BamMultiReader::MergeOrder, GetMergeOrder(), GetNextAlignment(), GetNextAlignmentCore() +*/ +void BamMultiReader::SetExplicitMergeOrder(BamMultiReader::MergeOrder order) { + d->SetExplicitMergeOrder(order); } -// updates the reference id stored in the BamMultiReader -// to reflect the current state of the readers -void BamMultiReader::UpdateReferenceID(void) { - // the alignments are sorted by position, so the first alignment will always have the lowest reference ID - if (alignments.begin()->second.second->RefID != CurrentRefID) { - // get the next reference id - // while there aren't any readers at the next ref id - // increment the ref id - int nextRefID = CurrentRefID; - while (alignments.begin()->second.second->RefID != nextRefID) { - ++nextRefID; - } - //cerr << "updating reference id from " << CurrentRefID << " to " << nextRefID << endl; - CurrentRefID = nextRefID; - } +/*! \fn bool BamMultiReader::SetRegion(const BamRegion& region) + \brief Sets a target region of interest + + Equivalent to calling BamReader::SetRegion() on all open BAM files. + + \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 ALL readers set the region successfully + \sa HasIndexes(), Jump(), BamReader::SetRegion() +*/ +bool BamMultiReader::SetRegion(const BamRegion& region) { + return d->SetRegion(region); } -// ValidateReaders checks that all the readers point to BAM files representing -// alignments against the same set of reference sequences, and that the -// sequences are identically ordered. If these checks fail the operation of -// the multireader is undefined, so we force program exit. -void BamMultiReader::ValidateReaders(void) const { - int firstRefCount = readers.front().first->GetReferenceCount(); - BamTools::RefVector firstRefData = readers.front().first->GetReferenceData(); - for (vector >::const_iterator it = readers.begin(); it != readers.end(); ++it) { - BamReader* reader = it->first; - BamTools::RefVector currentRefData = reader->GetReferenceData(); - BamTools::RefVector::const_iterator f = firstRefData.begin(); - BamTools::RefVector::const_iterator c = currentRefData.begin(); - if (reader->GetReferenceCount() != firstRefCount || firstRefData.size() != currentRefData.size()) { - cerr << "ERROR: mismatched number of references in " << reader->GetFilename() - << " expected " << firstRefCount - << " reference sequences but only found " << reader->GetReferenceCount() << endl; - exit(1); - } - // this will be ok; we just checked above that we have identically-sized sets of references - // here we simply check if they are all, in fact, equal in content - while (f != firstRefData.end()) { - if (f->RefName != c->RefName || f->RefLength != c->RefLength) { - cerr << "ERROR: mismatched references found in " << reader->GetFilename() - << " expected: " << endl; - for (BamTools::RefVector::const_iterator a = firstRefData.begin(); a != firstRefData.end(); ++a) - cerr << a->RefName << " " << a->RefLength << endl; - cerr << "but found: " << endl; - for (BamTools::RefVector::const_iterator a = currentRefData.begin(); a != currentRefData.end(); ++a) - cerr << a->RefName << " " << a->RefLength << endl; - exit(1); - } - ++f; ++c; - } - } +/*! \fn bool BamMultiReader::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. Equivalent to calling BamReader::SetRegion() on all open BAM files. + + \warning This function now 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 ALL readers set the region successfully + \sa HasIndexes(), Jump(), BamReader::SetRegion() +*/ +bool BamMultiReader::SetRegion(const int& leftRefID, + const int& leftPosition, + const int& rightRefID, + const int& rightPosition) +{ + return d->SetRegion( BamRegion(leftRefID, leftPosition, rightRefID, rightPosition) ); }