// ***************************************************************************
// BamMultiReader.cpp (c) 2010 Erik Garrison, Derek Barnett
// Marth Lab, Department of Biology, Boston College
-// All rights reserved.
// ---------------------------------------------------------------------------
-// Last modified: 18 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.
// ***************************************************************************
-#include <algorithm>
-#include <fstream>
-#include <iostream>
-#include <iterator>
-#include <sstream>
+#include "api/BamMultiReader.h"
+#include "api/internal/bam/BamMultiReader_p.h"
+using namespace BamTools;
+
#include <string>
#include <vector>
-#include "BGZF.h"
-#include "BamMultiReader.h"
-using namespace BamTools;
using namespace std;
-// -----------------------------------------------------
-// BamMultiReader implementation
-// -----------------------------------------------------
+/*! \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
+*/
+
+
-// 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<pair<BamReader*, BamAlignment*> >::iterator readerIter = readers.begin();
- vector<pair<BamReader*, BamAlignment*> >::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 useStandardIndex) {
- bool result = true;
- for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
- BamReader* reader = it->first;
- result &= reader->CreateIndex(useStandardIndex);
- }
- 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<string, bool> readGroups;
-
- // foreach extraction entry (each BAM file)
- for (vector<pair<BamReader*, BamAlignment*> >::const_iterator rs = readers.begin(); rs != readers.end(); ++rs) {
-
- BamReader* reader = rs->first;
- string headerText = reader->GetHeaderText();
- if ( headerText.empty() ) continue;
-
- map<string, bool> currentFileReadGroups;
- stringstream header(headerText);
- vector<string> lines;
- string item;
- while (getline(header, item))
- lines.push_back(item);
-
- for (vector<string>::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<std::string> 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<std::string> 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<pair<BamReader*, BamAlignment*> >::const_iterator readerIter = readers.begin();
- vector<pair<BamReader*, BamAlignment*> >::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<pair<BamReader*, BamAlignment*> >::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<string>& filenames, bool openIndexes, bool coreMode, bool useDefaultIndex) {
-
- // for filename in filenames
- fileNames = filenames; // save filenames in our multireader
- for (vector<string>::const_iterator it = filenames.begin(); it != filenames.end(); ++it) {
-
- const string filename = *it;
- BamReader* reader = new BamReader;
-
- bool openedOK = true;
- if (openIndexes) {
-
- // leave index filename empty
- // this allows BamReader & BamIndex to search for any available
- // useDefaultIndex gives hint to prefer BAI over BTI
- openedOK = reader->Open(filename, "", true, useDefaultIndex);
- }
-
- // 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 further error handling when openedOK is 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<pair<BamReader*, BamAlignment*> >::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<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
- BamReader* reader = it->first;
- result &= reader->Rewind();
- }
- return result;
+/*! \fn bool BamMultiReader::Open(const std::vector<std::string>& 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<std::string>& 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<std::string>& 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<std::string>& 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<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
- if (!it->first->SetRegion(region)) {
- cerr << "ERROR: could not jump " << it->first->GetFilename() << " to "
- << region.LeftRefID << ":" << region.LeftPosition
- << ".." << region.RightRefID << ":" << region.RightPosition << endl;
- }
- }
+ 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<pair<BamReader*, BamAlignment*> >::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<pair<BamReader*, BamAlignment*> >::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) );
}
projects / bamtools.git / blobdiff