1 // ***************************************************************************
2 // BamMultiReader_p.cpp (c) 2010 Derek Barnett, Erik Garrison
3 // Marth Lab, Department of Biology, Boston College
4 // ---------------------------------------------------------------------------
5 // Last modified: 14 January 2013 (DB)
6 // ---------------------------------------------------------------------------
7 // Functionality for simultaneously reading multiple BAM files
8 // *************************************************************************
10 #include "api/BamAlignment.h"
11 #include "api/BamMultiReader.h"
12 #include "api/SamConstants.h"
13 #include "api/algorithms/Sort.h"
14 #include "api/internal/bam/BamMultiReader_p.h"
15 using namespace BamTools;
16 using namespace BamTools::Internal;
26 BamMultiReaderPrivate::BamMultiReaderPrivate(void)
28 , m_hasUserMergeOrder(false)
29 , m_mergeOrder(BamMultiReader::RoundRobinMerge)
33 BamMultiReaderPrivate::~BamMultiReaderPrivate(void) {
37 // close all BAM files
38 bool BamMultiReaderPrivate::Close(void) {
40 m_errorString.clear();
42 if ( CloseFiles(Filenames()) )
45 const string currentError = m_errorString;
46 const string message = string("error encountered while closing all files: \n\t") + currentError;
47 SetErrorString("BamMultiReader::Close", message);
52 // close requested BAM file
53 bool BamMultiReaderPrivate::CloseFile(const string& filename) {
55 m_errorString.clear();
57 vector<string> filenames(1, filename);
58 if ( CloseFiles(filenames) )
61 const string currentError = m_errorString;
62 const string message = string("error while closing file: ") + filename + "\n" + currentError;
63 SetErrorString("BamMultiReader::CloseFile", message);
68 // close requested BAM files
69 bool BamMultiReaderPrivate::CloseFiles(const vector<string>& filenames) {
71 bool errorsEncountered = false;
72 m_errorString.clear();
74 // iterate over filenames
75 vector<string>::const_iterator filesIter = filenames.begin();
76 vector<string>::const_iterator filesEnd = filenames.end();
77 for ( ; filesIter != filesEnd; ++filesIter ) {
78 const string& filename = (*filesIter);
79 if ( filename.empty() ) continue;
81 // iterate over readers
82 vector<MergeItem>::iterator readerIter = m_readers.begin();
83 vector<MergeItem>::iterator readerEnd = m_readers.end();
84 for ( ; readerIter != readerEnd; ++readerIter ) {
85 MergeItem& item = (*readerIter);
86 BamReader* reader = item.Reader;
87 if ( reader == 0 ) continue;
89 // if reader matches requested filename
90 if ( reader->GetFilename() == filename ) {
92 // remove reader's entry from alignment cache
93 m_alignmentCache->Remove(reader);
95 // clean up reader & its alignment
96 if ( !reader->Close() ) {
97 m_errorString.append(1, '\t');
98 m_errorString.append(reader->GetErrorString());
99 m_errorString.append(1, '\n');
100 errorsEncountered = true;
105 // delete reader's alignment entry
106 BamAlignment* alignment = item.Alignment;
110 // remove reader from reader list
111 m_readers.erase(readerIter);
113 // on match, just go on to next filename
114 // (no need to keep looking and item iterator is invalid now anyway)
120 // make sure we clean up properly if all readers were closed
121 if ( m_readers.empty() ) {
124 if ( m_alignmentCache ) {
125 m_alignmentCache->Clear();
126 delete m_alignmentCache;
127 m_alignmentCache = 0;
131 m_hasUserMergeOrder = false;
132 m_mergeOrder = BamMultiReader::RoundRobinMerge;
135 // return whether all readers closed OK
136 return !errorsEncountered;
139 // creates index files for BAM files that don't have them
140 bool BamMultiReaderPrivate::CreateIndexes(const BamIndex::IndexType& type) {
142 bool errorsEncountered = false;
143 m_errorString.clear();
145 // iterate over readers
146 vector<MergeItem>::iterator itemIter = m_readers.begin();
147 vector<MergeItem>::iterator itemEnd = m_readers.end();
148 for ( ; itemIter != itemEnd; ++itemIter ) {
149 MergeItem& item = (*itemIter);
150 BamReader* reader = item.Reader;
151 if ( reader == 0 ) continue;
153 // if reader doesn't have an index, create one
154 if ( !reader->HasIndex() ) {
155 if ( !reader->CreateIndex(type) ) {
156 m_errorString.append(1, '\t');
157 m_errorString.append(reader->GetErrorString());
158 m_errorString.append(1, '\n');
159 errorsEncountered = true;
164 // check for errors encountered before returning success/fail
165 if ( errorsEncountered ) {
166 const string currentError = m_errorString;
167 const string message = string("error while creating index files: ") + "\n" + currentError;
168 SetErrorString("BamMultiReader::CreateIndexes", message);
174 IMultiMerger* BamMultiReaderPrivate::CreateAlignmentCache(void) {
176 // if no merge order set explicitly, use SAM header to lookup proper order
177 if ( !m_hasUserMergeOrder ) {
179 // fetch SamHeader from BAM files
180 SamHeader header = GetHeader();
182 // if BAM files are sorted by position
183 if ( header.SortOrder == Constants::SAM_HD_SORTORDER_COORDINATE )
184 m_mergeOrder = BamMultiReader::MergeByCoordinate;
186 // if BAM files are sorted by read name
187 if ( header.SortOrder == Constants::SAM_HD_SORTORDER_QUERYNAME )
188 m_mergeOrder = BamMultiReader::MergeByName;
190 // otherwise, sorting is either "unknown" or marked as "unsorted"
192 m_mergeOrder = BamMultiReader::RoundRobinMerge;
195 // use current merge order to create proper 'multi-merger'
196 switch ( m_mergeOrder ) {
198 // merge BAM files by position
199 case BamMultiReader::MergeByCoordinate :
200 return new MultiMerger<Algorithms::Sort::ByPosition>();
202 // merge BAM files by read name
203 case BamMultiReader::MergeByName :
204 return new MultiMerger<Algorithms::Sort::ByName>();
206 // sorting is "unknown", "unsorted" or "ignored"... so use unsorted merger
207 case BamMultiReader::RoundRobinMerge :
208 return new MultiMerger<Algorithms::Sort::Unsorted>();
210 // unknown merge order, can't create merger
216 const vector<string> BamMultiReaderPrivate::Filenames(void) const {
218 // init filename container
219 vector<string> filenames;
220 filenames.reserve( m_readers.size() );
222 // iterate over readers
223 vector<MergeItem>::const_iterator itemIter = m_readers.begin();
224 vector<MergeItem>::const_iterator itemEnd = m_readers.end();
225 for ( ; itemIter != itemEnd; ++itemIter ) {
226 const MergeItem& item = (*itemIter);
227 const BamReader* reader = item.Reader;
228 if ( reader == 0 ) continue;
230 // store filename if not empty
231 const string& filename = reader->GetFilename();
232 if ( !filename.empty() )
233 filenames.push_back(filename);
240 string BamMultiReaderPrivate::GetErrorString(void) const {
241 return m_errorString;
244 SamHeader BamMultiReaderPrivate::GetHeader(void) const {
245 const string& text = GetHeaderText();
246 return SamHeader(text);
249 // makes a virtual, unified header for all the bam files in the multireader
250 string BamMultiReaderPrivate::GetHeaderText(void) const {
252 // N.B. - right now, simply copies all header data from first BAM,
253 // and then appends RG's from other BAM files
254 // TODO: make this more intelligent wrt other header lines/fields
256 // if no readers open
257 const size_t numReaders = m_readers.size();
258 if ( numReaders == 0 ) return string();
260 // retrieve first reader's header
261 const MergeItem& firstItem = m_readers.front();
262 const BamReader* reader = firstItem.Reader;
263 if ( reader == 0 ) return string();
264 SamHeader mergedHeader = reader->GetHeader();
266 // iterate over any remaining readers (skipping the first)
267 for ( size_t i = 1; i < numReaders; ++i ) {
268 const MergeItem& item = m_readers.at(i);
269 const BamReader* reader = item.Reader;
270 if ( reader == 0 ) continue;
272 // retrieve current reader's header
273 const SamHeader currentHeader = reader->GetHeader();
275 // append current reader's RG entries to merged header
276 // N.B. - SamReadGroupDictionary handles duplicate-checking
277 mergedHeader.ReadGroups.Add(currentHeader.ReadGroups);
279 // TODO: merge anything else??
282 // return stringified header
283 return mergedHeader.ToString();
286 BamMultiReader::MergeOrder BamMultiReaderPrivate::GetMergeOrder(void) const {
290 // get next alignment among all files
291 bool BamMultiReaderPrivate::GetNextAlignment(BamAlignment& al) {
292 return PopNextCachedAlignment(al, true);
295 // get next alignment among all files without parsing character data from alignments
296 bool BamMultiReaderPrivate::GetNextAlignmentCore(BamAlignment& al) {
297 return PopNextCachedAlignment(al, false);
300 // ---------------------------------------------------------------------------------------
302 // NB: The following GetReferenceX() functions assume that we have identical
303 // references for all BAM files. We enforce this by invoking the
304 // ValidateReaders() method to verify that our reference data is the same
305 // across all files on Open - so we will not encounter a situation in which
306 // there is a mismatch and we are still live.
308 // ---------------------------------------------------------------------------------------
310 // returns the number of reference sequences
311 int BamMultiReaderPrivate::GetReferenceCount(void) const {
313 // handle empty multireader
314 if ( m_readers.empty() ) return 0;
316 // return reference count from first reader
317 const MergeItem& item = m_readers.front();
318 const BamReader* reader = item.Reader;
319 if ( reader == 0 ) return 0;
321 return reader->GetReferenceCount();
324 // returns vector of reference objects
325 const RefVector BamMultiReaderPrivate::GetReferenceData(void) const {
327 // handle empty multireader
328 if ( m_readers.empty() ) return RefVector();
330 // return reference data from first BamReader
331 const MergeItem& item = m_readers.front();
332 const BamReader* reader = item.Reader;
333 if ( reader == 0 ) return RefVector();
335 return reader->GetReferenceData();
338 // returns refID from reference name
339 int BamMultiReaderPrivate::GetReferenceID(const string& refName) const {
341 // handle empty multireader
342 if ( m_readers.empty() ) return -1;
344 // return reference ID from first BamReader
345 const MergeItem& item = m_readers.front();
346 const BamReader* reader = item.Reader;
347 if ( reader == 0 ) return -1;
349 return reader->GetReferenceID(refName);
351 // ---------------------------------------------------------------------------------------
353 // returns true if all readers have index data available
354 // this is useful to indicate whether Jump() or SetRegion() are possible
355 bool BamMultiReaderPrivate::HasIndexes(void) const {
357 // handle empty multireader
358 if ( m_readers.empty() )
363 // iterate over readers
364 vector<MergeItem>::const_iterator readerIter = m_readers.begin();
365 vector<MergeItem>::const_iterator readerEnd = m_readers.end();
366 for ( ; readerIter != readerEnd; ++readerIter ) {
367 const MergeItem& item = (*readerIter);
368 const BamReader* reader = item.Reader;
369 if ( reader == 0 ) continue;
371 // see if current reader has index data
372 result &= reader->HasIndex();
378 // returns true if multireader has open readers
379 bool BamMultiReaderPrivate::HasOpenReaders(void) {
381 // iterate over readers
382 vector<MergeItem>::const_iterator readerIter = m_readers.begin();
383 vector<MergeItem>::const_iterator readerEnd = m_readers.end();
384 for ( ; readerIter != readerEnd; ++readerIter ) {
385 const MergeItem& item = (*readerIter);
386 const BamReader* reader = item.Reader;
387 if ( reader == 0 ) continue;
389 // return true whenever an open reader is found
390 if ( reader->IsOpen() ) return true;
397 // performs random-access jump using (refID, position) as a left-bound
398 bool BamMultiReaderPrivate::Jump(int refID, int position) {
400 // NB: While it may make sense to track readers in which we can
401 // successfully Jump, in practice a failure of Jump means "no
402 // alignments here." It makes sense to simply accept the failure,
403 // UpdateAlignments(), and continue.
405 // iterate over readers
406 vector<MergeItem>::iterator readerIter = m_readers.begin();
407 vector<MergeItem>::iterator readerEnd = m_readers.end();
408 for ( ; readerIter != readerEnd; ++readerIter ) {
409 MergeItem& item = (*readerIter);
410 BamReader* reader = item.Reader;
411 if ( reader == 0 ) continue;
413 // jump in each BamReader to position of interest
414 reader->Jump(refID, position);
417 // returns status of cache update
418 return UpdateAlignmentCache();
421 // locate (& load) index files for BAM readers that don't already have one loaded
422 bool BamMultiReaderPrivate::LocateIndexes(const BamIndex::IndexType& preferredType) {
424 bool errorsEncountered = false;
425 m_errorString.clear();
427 // iterate over readers
428 vector<MergeItem>::iterator readerIter = m_readers.begin();
429 vector<MergeItem>::iterator readerEnd = m_readers.end();
430 for ( ; readerIter != readerEnd; ++readerIter ) {
431 MergeItem& item = (*readerIter);
432 BamReader* reader = item.Reader;
433 if ( reader == 0 ) continue;
435 // if reader has no index, try to locate one
436 if ( !reader->HasIndex() ) {
437 if ( !reader->LocateIndex(preferredType) ) {
438 m_errorString.append(1, '\t');
439 m_errorString.append(reader->GetErrorString());
440 m_errorString.append(1, '\n');
441 errorsEncountered = true;
446 // check for errors encountered before returning success/fail
447 if ( errorsEncountered ) {
448 const string currentError = m_errorString;
449 const string message = string("error while locating index files: ") + "\n" + currentError;
450 SetErrorString("BamMultiReader::LocatingIndexes", message);
457 bool BamMultiReaderPrivate::Open(const vector<string>& filenames) {
459 m_errorString.clear();
461 // put all current readers back at beginning (refreshes alignment cache)
463 const string currentError = m_errorString;
464 const string message = string("unable to rewind existing readers: \n\t") + currentError;
465 SetErrorString("BamMultiReader::Open", message);
469 // iterate over filenames
470 bool errorsEncountered = false;
471 vector<string>::const_iterator filenameIter = filenames.begin();
472 vector<string>::const_iterator filenameEnd = filenames.end();
473 for ( ; filenameIter != filenameEnd; ++filenameIter ) {
474 const string& filename = (*filenameIter);
475 if ( filename.empty() ) continue;
477 // attempt to open BamReader
478 BamReader* reader = new BamReader;
479 const bool readerOpened = reader->Open(filename);
481 // if opened OK, store it
483 m_readers.push_back( MergeItem(reader, new BamAlignment) );
485 // otherwise store error & clean up invalid reader
487 m_errorString.append(1, '\t');
488 m_errorString += string("unable to open file: ") + filename;
489 m_errorString.append(1, '\n');
490 errorsEncountered = true;
497 // check for errors while opening
498 if ( errorsEncountered ) {
499 const string currentError = m_errorString;
500 const string message = string("unable to open all files: \t\n") + currentError;
501 SetErrorString("BamMultiReader::Open", message);
505 // check for BAM file consistency
506 if ( !ValidateReaders() ) {
507 const string currentError = m_errorString;
508 const string message = string("unable to open inconsistent files: \t\n") + currentError;
509 SetErrorString("BamMultiReader::Open", message);
513 // update alignment cache
514 return UpdateAlignmentCache();
517 bool BamMultiReaderPrivate::OpenFile(const std::string& filename) {
518 vector<string> filenames(1, filename);
519 if ( Open(filenames) )
522 const string currentError = m_errorString;
523 const string message = string("could not open file: ") + filename + "\n\t" + currentError;
524 SetErrorString("BamMultiReader::OpenFile", message);
529 bool BamMultiReaderPrivate::OpenIndexes(const vector<string>& indexFilenames) {
531 // TODO: This needs to be cleaner - should not assume same order.
532 // And either way, shouldn't start at first reader. Should start at
533 // first reader without an index?
535 // make sure same number of index filenames as readers
536 if ( m_readers.size() != indexFilenames.size() ) {
537 const string message("size of index file list does not match current BAM file count");
538 SetErrorString("BamMultiReader::OpenIndexes", message);
542 bool errorsEncountered = false;
543 m_errorString.clear();
545 // iterate over BamReaders
546 vector<string>::const_iterator indexFilenameIter = indexFilenames.begin();
547 vector<string>::const_iterator indexFilenameEnd = indexFilenames.end();
548 vector<MergeItem>::iterator readerIter = m_readers.begin();
549 vector<MergeItem>::iterator readerEnd = m_readers.end();
550 for ( ; readerIter != readerEnd; ++readerIter ) {
551 MergeItem& item = (*readerIter);
552 BamReader* reader = item.Reader;
554 // open index filename on reader
556 const string& indexFilename = (*indexFilenameIter);
557 if ( !reader->OpenIndex(indexFilename) ) {
558 m_errorString.append(1, '\t');
559 m_errorString += reader->GetErrorString();
560 m_errorString.append(1, '\n');
561 errorsEncountered = true;
565 // increment filename iterator, skip if no more index files to open
566 if ( ++indexFilenameIter == indexFilenameEnd )
570 // return success/fail
571 if ( errorsEncountered ) {
572 const string currentError = m_errorString;
573 const string message = string("could not open all index files: \n\t") + currentError;
574 SetErrorString("BamMultiReader::OpenIndexes", message);
580 bool BamMultiReaderPrivate::PopNextCachedAlignment(BamAlignment& al, const bool needCharData) {
582 // skip if no alignments available
583 if ( m_alignmentCache == 0 || m_alignmentCache->IsEmpty() )
586 // pop next merge item entry from cache
587 MergeItem item = m_alignmentCache->TakeFirst();
588 BamReader* reader = item.Reader;
589 BamAlignment* alignment = item.Alignment;
590 if ( reader == 0 || alignment == 0 )
593 // set char data if requested
594 if ( needCharData ) {
595 alignment->BuildCharData();
596 alignment->Filename = reader->GetFilename();
599 // store cached alignment into destination parameter (by copy)
602 // load next alignment from reader & store in cache
603 SaveNextAlignment(reader, alignment);
607 // returns BAM file pointers to beginning of alignment data & resets alignment cache
608 bool BamMultiReaderPrivate::Rewind(void) {
610 // skip if no readers open
611 if ( m_readers.empty() )
614 // attempt to rewind files
615 if ( !RewindReaders() ) {
616 const string currentError = m_errorString;
617 const string message = string("could not rewind readers: \n\t") + currentError;
618 SetErrorString("BamMultiReader::Rewind", message);
622 // return status of cache update
623 return UpdateAlignmentCache();
626 // returns BAM file pointers to beginning of alignment data
627 bool BamMultiReaderPrivate::RewindReaders(void) {
629 m_errorString.clear();
630 bool errorsEncountered = false;
632 // iterate over readers
633 vector<MergeItem>::iterator readerIter = m_readers.begin();
634 vector<MergeItem>::iterator readerEnd = m_readers.end();
635 for ( ; readerIter != readerEnd; ++readerIter ) {
636 MergeItem& item = (*readerIter);
637 BamReader* reader = item.Reader;
638 if ( reader == 0 ) continue;
640 // attempt rewind on BamReader
641 if ( !reader->Rewind() ) {
642 m_errorString.append(1, '\t');
643 m_errorString.append( reader->GetErrorString() );
644 m_errorString.append(1, '\n');
645 errorsEncountered = true;
649 return !errorsEncountered;
652 void BamMultiReaderPrivate::SaveNextAlignment(BamReader* reader, BamAlignment* alignment) {
654 // if can read alignment from reader, store in cache
656 // N.B. - lazy building of alignment's char data - populated only:
657 // automatically by alignment cache to maintain its sorting OR
658 // on demand from client call to future call to GetNextAlignment()
660 if ( reader->GetNextAlignmentCore(*alignment) )
661 m_alignmentCache->Add( MergeItem(reader, alignment) );
664 bool BamMultiReaderPrivate::SetExplicitMergeOrder(BamMultiReader::MergeOrder order) {
666 // set new merge flags
667 m_hasUserMergeOrder = true;
668 m_mergeOrder = order;
670 // remove any existing merger (storing any existing data sitting in the cache)
671 vector<MergeItem> currentCacheData;
672 if ( m_alignmentCache ) {
673 while ( !m_alignmentCache->IsEmpty() )
674 currentCacheData.push_back( m_alignmentCache->TakeFirst() );
675 delete m_alignmentCache;
676 m_alignmentCache = 0;
679 // create new cache using the new merge flags
680 m_alignmentCache = CreateAlignmentCache();
681 if ( m_alignmentCache == 0 ) {
682 SetErrorString("BamMultiReader::SetExplicitMergeOrder", "requested order is unrecognized");
686 // push current data onto new cache
687 vector<MergeItem>::const_iterator readerIter = currentCacheData.begin();
688 vector<MergeItem>::const_iterator readerEnd = currentCacheData.end();
689 for ( ; readerIter != readerEnd; ++readerIter ) {
690 const MergeItem& item = (*readerIter);
691 m_alignmentCache->Add(item);
698 void BamMultiReaderPrivate::SetErrorString(const string& where, const string& what) const {
699 static const string SEPARATOR = ": ";
700 m_errorString = where + SEPARATOR + what;
703 bool BamMultiReaderPrivate::SetRegion(const BamRegion& region) {
705 // NB: While it may make sense to track readers in which we can
706 // successfully SetRegion, In practice a failure of SetRegion means "no
707 // alignments here." It makes sense to simply accept the failure,
708 // UpdateAlignments(), and continue.
710 // iterate over alignments
711 vector<MergeItem>::iterator readerIter = m_readers.begin();
712 vector<MergeItem>::iterator readerEnd = m_readers.end();
713 for ( ; readerIter != readerEnd; ++readerIter ) {
714 MergeItem& item = (*readerIter);
715 BamReader* reader = item.Reader;
716 if ( reader == 0 ) continue;
718 // set region of interest
719 reader->SetRegion(region);
722 // return status of cache update
723 return UpdateAlignmentCache();
726 // updates our alignment cache
727 bool BamMultiReaderPrivate::UpdateAlignmentCache(void) {
729 // create alignment cache if not created yet
730 if ( m_alignmentCache == 0 ) {
731 m_alignmentCache = CreateAlignmentCache();
732 if ( m_alignmentCache == 0 ) {
733 SetErrorString("BamMultiReader::UpdateAlignmentCache", "unable to create new alignment cache");
738 // clear any prior cache data
739 m_alignmentCache->Clear();
741 // iterate over readers
742 vector<MergeItem>::iterator readerIter = m_readers.begin();
743 vector<MergeItem>::iterator readerEnd = m_readers.end();
744 for ( ; readerIter != readerEnd; ++readerIter ) {
745 MergeItem& item = (*readerIter);
746 BamReader* reader = item.Reader;
747 BamAlignment* alignment = item.Alignment;
748 if ( reader == 0 || alignment == 0 ) continue;
750 // save next alignment from each reader in cache
751 SaveNextAlignment(reader, alignment);
754 // if we get here, ok
758 // ValidateReaders checks that all the readers point to BAM files representing
759 // alignments against the same set of reference sequences, and that the
760 // sequences are identically ordered. If these checks fail the operation of
761 // the multireader is undefined, so we force program exit.
762 bool BamMultiReaderPrivate::ValidateReaders(void) const {
764 m_errorString.clear();
766 // skip if 0 or 1 readers opened
767 if ( m_readers.empty() || (m_readers.size() == 1) )
770 // retrieve first reader
771 const MergeItem& firstItem = m_readers.front();
772 const BamReader* firstReader = firstItem.Reader;
773 if ( firstReader == 0 ) return false;
775 // retrieve first reader's header data
776 const SamHeader& firstReaderHeader = firstReader->GetHeader();
777 const string& firstReaderSortOrder = firstReaderHeader.SortOrder;
779 // retrieve first reader's reference data
780 const RefVector& firstReaderRefData = firstReader->GetReferenceData();
781 const int firstReaderRefCount = firstReader->GetReferenceCount();
782 const int firstReaderRefSize = firstReaderRefData.size();
784 // iterate over all readers
785 vector<MergeItem>::const_iterator readerIter = m_readers.begin();
786 vector<MergeItem>::const_iterator readerEnd = m_readers.end();
787 for ( ; readerIter != readerEnd; ++readerIter ) {
788 const MergeItem& item = (*readerIter);
789 BamReader* reader = item.Reader;
790 if ( reader == 0 ) continue;
792 // get current reader's header data
793 const SamHeader& currentReaderHeader = reader->GetHeader();
794 const string& currentReaderSortOrder = currentReaderHeader.SortOrder;
796 // check compatible sort order
797 if ( currentReaderSortOrder != firstReaderSortOrder ) {
798 const string message = string("mismatched sort order in ") + reader->GetFilename() +
799 ", expected " + firstReaderSortOrder +
800 ", but found " + currentReaderSortOrder;
801 SetErrorString("BamMultiReader::ValidateReaders", message);
805 // get current reader's reference data
806 const RefVector currentReaderRefData = reader->GetReferenceData();
807 const int currentReaderRefCount = reader->GetReferenceCount();
808 const int currentReaderRefSize = currentReaderRefData.size();
810 // init reference data iterators
811 RefVector::const_iterator firstRefIter = firstReaderRefData.begin();
812 RefVector::const_iterator firstRefEnd = firstReaderRefData.end();
813 RefVector::const_iterator currentRefIter = currentReaderRefData.begin();
815 // compare reference counts from BamReader ( & container size, in case of BR error)
816 if ( (currentReaderRefCount != firstReaderRefCount) ||
817 (firstReaderRefSize != currentReaderRefSize) )
820 s << "mismatched reference count in " << reader->GetFilename()
821 << ", expected " << firstReaderRefCount
822 << ", but found " << currentReaderRefCount;
823 SetErrorString("BamMultiReader::ValidateReaders", s.str());
827 // this will be ok; we just checked above that we have identically-sized sets of references
828 // here we simply check if they are all, in fact, equal in content
829 while ( firstRefIter != firstRefEnd ) {
830 const RefData& firstRef = (*firstRefIter);
831 const RefData& currentRef = (*currentRefIter);
833 // compare reference name & length
834 if ( (firstRef.RefName != currentRef.RefName) ||
835 (firstRef.RefLength != currentRef.RefLength) )
838 s << "mismatched references found in" << reader->GetFilename()
839 << "expected: " << endl;
841 // print first reader's reference data
842 RefVector::const_iterator refIter = firstReaderRefData.begin();
843 RefVector::const_iterator refEnd = firstReaderRefData.end();
844 for ( ; refIter != refEnd; ++refIter ) {
845 const RefData& entry = (*refIter);
847 s << entry.RefName << " " << endl;
850 s << "but found: " << endl;
852 // print current reader's reference data
853 refIter = currentReaderRefData.begin();
854 refEnd = currentReaderRefData.end();
855 for ( ; refIter != refEnd; ++refIter ) {
856 const RefData& entry = (*refIter);
857 s << entry.RefName << " " << entry.RefLength << endl;
860 SetErrorString("BamMultiReader::ValidateReaders", s.str());
870 // if we get here, everything checks out