src/api/internal/BamMultiReader_p.cpp

   1 // ***************************************************************************
   2 // BamMultiReader_p.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 3 October 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Functionality for simultaneously reading multiple BAM files
   8 // *************************************************************************
   9
  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/BamMultiReader_p.h>
  15 using namespace BamTools;
  16 using namespace BamTools::Internal;
  17
  18 #include <algorithm>
  19 #include <fstream>
  20 #include <iostream>
  21 #include <iterator>
  22 #include <sstream>
  23 using namespace std;
  24
  25 // ctor
  26 BamMultiReaderPrivate::BamMultiReaderPrivate(void)
  27     : m_alignmentCache(0)
  28 { }
  29
  30 // dtor
  31 BamMultiReaderPrivate::~BamMultiReaderPrivate(void) {
  32
  33     // close all open BAM readers (& clean up cache)
  34     Close();
  35 }
  36
  37 // close all BAM files
  38 void BamMultiReaderPrivate::Close(void) {
  39     CloseFiles( Filenames() );
  40 }
  41
  42 // close requested BAM file
  43 void BamMultiReaderPrivate::CloseFile(const string& filename) {
  44     vector<string> filenames(1, filename);
  45     CloseFiles(filenames);
  46 }
  47
  48 // close requested BAM files
  49 void BamMultiReaderPrivate::CloseFiles(const vector<string>& filenames) {
  50
  51     // iterate over filenames
  52     vector<string>::const_iterator filesIter = filenames.begin();
  53     vector<string>::const_iterator filesEnd  = filenames.end();
  54     for ( ; filesIter != filesEnd; ++filesIter ) {
  55         const string& filename = (*filesIter);
  56         if ( filename.empty() ) continue;
  57
  58         // iterate over readers
  59         vector<MergeItem>::iterator readerIter = m_readers.begin();
  60         vector<MergeItem>::iterator readerEnd  = m_readers.end();
  61         for ( ; readerIter != readerEnd; ++readerIter ) {
  62             MergeItem& item = (*readerIter);
  63             BamReader* reader = item.Reader;
  64             if ( reader == 0 ) continue;
  65
  66             // if reader matches requested filename
  67             if ( reader->GetFilename() == filename ) {
  68
  69                 // remove reader's entry from alignment cache
  70                 m_alignmentCache->Remove(reader);
  71
  72                 // clean up reader & its alignment
  73                 reader->Close();
  74                 delete reader;
  75                 reader = 0;
  76
  77                 // delete reader's alignment entry
  78                 BamAlignment* alignment = item.Alignment;
  79                 delete alignment;
  80                 alignment = 0;
  81
  82                 // remove reader from reader list
  83                 m_readers.erase(readerIter);
  84
  85                 // on match, just go on to next filename
  86                 // (no need to keep looking and item iterator is invalid now anyway)
  87                 break;
  88             }
  89         }
  90     }
  91
  92     // make sure alignment cache is cleaned up if all readers closed
  93     if ( m_readers.empty() && m_alignmentCache ) {
  94         m_alignmentCache->Clear();
  95         delete m_alignmentCache;
  96         m_alignmentCache = 0;
  97     }
  98 }
  99
 100 // creates index files for BAM files that don't have them
 101 bool BamMultiReaderPrivate::CreateIndexes(const BamIndex::IndexType& type) {
 102
 103     bool result = true;
 104
 105     // iterate over readers
 106     vector<MergeItem>::iterator itemIter = m_readers.begin();
 107     vector<MergeItem>::iterator itemEnd  = m_readers.end();
 108     for ( ; itemIter != itemEnd; ++itemIter ) {
 109         MergeItem& item = (*itemIter);
 110         BamReader* reader = item.Reader;
 111         if ( reader == 0 ) continue;
 112
 113         // if reader doesn't have an index, create one
 114         if ( !reader->HasIndex() )
 115             result &= reader->CreateIndex(type);
 116     }
 117
 118     return result;
 119 }
 120
 121 IMultiMerger* BamMultiReaderPrivate::CreateAlignmentCache(void) const {
 122
 123     // fetch SamHeader
 124     SamHeader header = GetHeader();
 125
 126     // if BAM files are sorted by position
 127     if ( header.SortOrder == Constants::SAM_HD_SORTORDER_COORDINATE )
 128         return new MultiMerger<Algorithms::Sort::ByPosition>();
 129
 130     // if BAM files are sorted by read name
 131     if ( header.SortOrder == Constants::SAM_HD_SORTORDER_QUERYNAME )
 132         return new MultiMerger<Algorithms::Sort::ByName>();
 133
 134     // otherwise "unknown" or "unsorted", use unsorted merger and just read in
 135     return new MultiMerger<Algorithms::Sort::Unsorted>();
 136 }
 137
 138 const vector<string> BamMultiReaderPrivate::Filenames(void) const {
 139
 140     // init filename container
 141     vector<string> filenames;
 142     filenames.reserve( m_readers.size() );
 143
 144     // iterate over readers
 145     vector<MergeItem>::const_iterator itemIter = m_readers.begin();
 146     vector<MergeItem>::const_iterator itemEnd  = m_readers.end();
 147     for ( ; itemIter != itemEnd; ++itemIter ) {
 148         const MergeItem& item = (*itemIter);
 149         const BamReader* reader = item.Reader;
 150         if ( reader == 0 ) continue;
 151
 152         // store filename if not empty
 153         const string& filename = reader->GetFilename();
 154         if ( !filename.empty() )
 155             filenames.push_back(filename);
 156     }
 157
 158     // return result
 159     return filenames;
 160 }
 161
 162 SamHeader BamMultiReaderPrivate::GetHeader(void) const {
 163     const string& text = GetHeaderText();
 164     return SamHeader(text);
 165 }
 166
 167 // makes a virtual, unified header for all the bam files in the multireader
 168 string BamMultiReaderPrivate::GetHeaderText(void) const {
 169
 170     // N.B. - right now, simply copies all header data from first BAM,
 171     //        and then appends RG's from other BAM files
 172     // TODO: make this more intelligent wrt other header lines/fields
 173
 174     // if no readers open
 175     const size_t numReaders = m_readers.size();
 176     if ( numReaders == 0 ) return string();
 177
 178     // retrieve first reader's header
 179     const MergeItem& firstItem = m_readers.front();
 180     const BamReader* reader = firstItem.Reader;
 181     if ( reader == 0 ) return string();
 182     SamHeader mergedHeader = reader->GetHeader();
 183
 184     // iterate over any remaining readers (skipping the first)
 185     for ( size_t i = 1; i < numReaders; ++i ) {
 186         const MergeItem& item = m_readers.at(i);
 187         const BamReader* reader = item.Reader;
 188         if ( reader == 0 ) continue;
 189
 190         // retrieve current reader's header
 191         const SamHeader currentHeader = reader->GetHeader();
 192
 193         // append current reader's RG entries to merged header
 194         // N.B. - SamReadGroupDictionary handles duplicate-checking
 195         mergedHeader.ReadGroups.Add(currentHeader.ReadGroups);
 196
 197         // TODO: merge anything else??
 198     }
 199
 200     // return stringified header
 201     return mergedHeader.ToString();
 202 }
 203
 204 // get next alignment among all files
 205 bool BamMultiReaderPrivate::GetNextAlignment(BamAlignment& al) {
 206     return PopNextCachedAlignment(al, true);
 207 }
 208
 209 // get next alignment among all files without parsing character data from alignments
 210 bool BamMultiReaderPrivate::GetNextAlignmentCore(BamAlignment& al) {
 211     return PopNextCachedAlignment(al, false);
 212 }
 213
 214 // ---------------------------------------------------------------------------------------
 215 //
 216 // NB: The following GetReferenceX() functions assume that we have identical
 217 // references for all BAM files.  We enforce this by invoking the
 218 // ValidateReaders() method to verify that our reference data is the same
 219 // across all files on Open - so we will not encounter a situation in which
 220 // there is a mismatch and we are still live.
 221 //
 222 // ---------------------------------------------------------------------------------------
 223
 224 // returns the number of reference sequences
 225 int BamMultiReaderPrivate::GetReferenceCount(void) const {
 226
 227     // handle empty multireader
 228     if ( m_readers.empty() )
 229         return 0;
 230
 231     // return reference count from first reader
 232     const MergeItem& item = m_readers.front();
 233     const BamReader* reader = item.Reader;
 234     if ( reader ) return reader->GetReferenceCount();
 235
 236     // invalid reader
 237     return 0;
 238 }
 239
 240 // returns vector of reference objects
 241 const RefVector BamMultiReaderPrivate::GetReferenceData(void) const {
 242
 243     // handle empty multireader
 244     if ( m_readers.empty() )
 245         return RefVector();
 246
 247     // return reference data from first BamReader
 248     const MergeItem& item = m_readers.front();
 249     const BamReader* reader = item.Reader;
 250     if ( reader ) return reader->GetReferenceData();
 251
 252     // invalid reader
 253     return RefVector();
 254 }
 255
 256 // returns refID from reference name
 257 int BamMultiReaderPrivate::GetReferenceID(const string& refName) const {
 258
 259     // handle empty multireader
 260     if ( m_readers.empty() )
 261         return -1;
 262
 263     // return reference ID from first BamReader
 264     const MergeItem& item = m_readers.front();
 265     const BamReader* reader = item.Reader;
 266     if ( reader ) return reader->GetReferenceID(refName);
 267
 268     // invalid reader
 269     return -1;
 270 }
 271 // ---------------------------------------------------------------------------------------
 272
 273 // returns true if all readers have index data available
 274 // this is useful to indicate whether Jump() or SetRegion() are possible
 275 bool BamMultiReaderPrivate::HasIndexes(void) const {
 276
 277     // handle empty multireader
 278     if ( m_readers.empty() )
 279         return false;
 280
 281     bool result = true;
 282
 283     // iterate over readers
 284     vector<MergeItem>::const_iterator readerIter = m_readers.begin();
 285     vector<MergeItem>::const_iterator readerEnd  = m_readers.end();
 286     for ( ; readerIter != readerEnd; ++readerIter ) {
 287         const MergeItem& item = (*readerIter);
 288         const BamReader* reader = item.Reader;
 289         if ( reader  == 0 ) continue;
 290
 291         // see if current reader has index data
 292         result &= reader->HasIndex();
 293     }
 294
 295     return result;
 296 }
 297
 298 // returns true if multireader has open readers
 299 bool BamMultiReaderPrivate::HasOpenReaders(void) {
 300
 301     // iterate over readers
 302     vector<MergeItem>::const_iterator readerIter = m_readers.begin();
 303     vector<MergeItem>::const_iterator readerEnd  = m_readers.end();
 304     for ( ; readerIter != readerEnd; ++readerIter ) {
 305         const MergeItem& item = (*readerIter);
 306         const BamReader* reader = item.Reader;
 307         if ( reader == 0 ) continue;
 308
 309         // return true whenever an open reader is found
 310         if ( reader->IsOpen() ) return true;
 311     }
 312
 313     // no readers open
 314     return false;
 315 }
 316
 317 // performs random-access jump using (refID, position) as a left-bound
 318 bool BamMultiReaderPrivate::Jump(int refID, int position) {
 319
 320     // NB: While it may make sense to track readers in which we can
 321     // successfully Jump, in practice a failure of Jump means "no
 322     // alignments here."  It makes sense to simply accept the failure,
 323     // UpdateAlignments(), and continue.
 324
 325     // iterate over readers
 326     vector<MergeItem>::iterator readerIter = m_readers.begin();
 327     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 328     for ( ; readerIter != readerEnd; ++readerIter ) {
 329         MergeItem& item = (*readerIter);
 330         BamReader* reader = item.Reader;
 331         if ( reader == 0 ) continue;
 332
 333         // attempt jump() on each
 334         if ( !reader->Jump(refID, position) ) {
 335             cerr << "BamMultiReader ERROR: could not jump " << reader->GetFilename()
 336                  << " to " << refID << ":" << position << endl;
 337         }
 338     }
 339
 340     // returns status of cache update
 341     return UpdateAlignmentCache();
 342 }
 343
 344 // locate (& load) index files for BAM readers that don't already have one loaded
 345 bool BamMultiReaderPrivate::LocateIndexes(const BamIndex::IndexType& preferredType) {
 346
 347     bool result = true;
 348
 349     // iterate over readers
 350     vector<MergeItem>::iterator readerIter = m_readers.begin();
 351     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 352     for ( ; readerIter != readerEnd; ++readerIter ) {
 353         MergeItem& item = (*readerIter);
 354         BamReader* reader = item.Reader;
 355         if ( reader == 0 ) continue;
 356
 357         // if reader has no index, try to locate one
 358         if ( !reader->HasIndex() )
 359             result &= reader->LocateIndex(preferredType);
 360     }
 361
 362     return result;
 363 }
 364
 365 // opens BAM files
 366 bool BamMultiReaderPrivate::Open(const vector<string>& filenames) {
 367
 368     bool openedOk = true;
 369
 370     // put all current readers back at beginning
 371     openedOk &= Rewind();
 372
 373     // iterate over filenames
 374     vector<string>::const_iterator filenameIter = filenames.begin();
 375     vector<string>::const_iterator filenameEnd  = filenames.end();
 376     for ( ; filenameIter != filenameEnd; ++filenameIter ) {
 377         const string& filename = (*filenameIter);
 378         if ( filename.empty() ) continue;
 379
 380         // attempt to open BamReader
 381         BamReader* reader = new BamReader;
 382         const bool readerOpened = reader->Open(filename);
 383
 384         // if opened OK, store it
 385         if ( readerOpened )
 386             m_readers.push_back( MergeItem(reader, new BamAlignment) );
 387
 388         // otherwise clean up invalid reader
 389         else delete reader;
 390
 391         // update method return status
 392         openedOk &= readerOpened;
 393     }
 394
 395     // if more than one reader open, check for consistency
 396     if ( m_readers.size() > 1 )
 397         openedOk &= ValidateReaders();
 398
 399     // update alignment cache
 400     openedOk &= UpdateAlignmentCache();
 401
 402     // return success
 403     return openedOk;
 404 }
 405
 406 bool BamMultiReaderPrivate::OpenFile(const std::string& filename) {
 407     vector<string> filenames(1, filename);
 408     return Open(filenames);
 409 }
 410
 411 bool BamMultiReaderPrivate::OpenIndexes(const vector<string>& indexFilenames) {
 412
 413     // TODO: This needs to be cleaner - should not assume same order.
 414     //       And either way, shouldn't start at first reader.  Should start at
 415     //       first reader without an index?
 416
 417     // make sure same number of index filenames as readers
 418     if ( m_readers.size() != indexFilenames.size() )
 419         return false;
 420
 421     // init result flag
 422     bool result = true;
 423
 424     // iterate over BamReaders
 425     vector<string>::const_iterator indexFilenameIter = indexFilenames.begin();
 426     vector<string>::const_iterator indexFilenameEnd  = indexFilenames.end();
 427     vector<MergeItem>::iterator readerIter = m_readers.begin();
 428     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 429     for ( ; readerIter != readerEnd; ++readerIter ) {
 430         MergeItem& item = (*readerIter);
 431         BamReader* reader = item.Reader;
 432
 433         // open index filename on reader
 434         if ( reader ) {
 435             const string& indexFilename = (*indexFilenameIter);
 436             result &= reader->OpenIndex(indexFilename);
 437         }
 438
 439         // increment filename iterator, skip if no more index files to open
 440         if ( ++indexFilenameIter == indexFilenameEnd )
 441             break;
 442     }
 443
 444     // TODO: any validation needed here??
 445
 446     // return success/fail
 447     return result;
 448 }
 449
 450 bool BamMultiReaderPrivate::PopNextCachedAlignment(BamAlignment& al, const bool needCharData) {
 451
 452     // skip if no alignments available
 453     if ( m_alignmentCache == 0 || m_alignmentCache->IsEmpty() )
 454         return false;
 455
 456     // pop next merge item entry from cache
 457     MergeItem item = m_alignmentCache->TakeFirst();
 458     BamReader* reader = item.Reader;
 459     BamAlignment* alignment = item.Alignment;
 460     if ( reader == 0 || alignment == 0 )
 461         return false;
 462
 463     // set char data if requested
 464     if ( needCharData ) {
 465         alignment->BuildCharData();
 466         alignment->Filename = reader->GetFilename();
 467     }
 468
 469     // store cached alignment into destination parameter (by copy)
 470     al = *alignment;
 471
 472     // load next alignment from reader & store in cache
 473     SaveNextAlignment(reader, alignment);
 474
 475     // return success
 476     return true;
 477 }
 478
 479 // returns BAM file pointers to beginning of alignment data & resets alignment cache
 480 bool BamMultiReaderPrivate::Rewind(void) {
 481
 482     // attempt to rewind files
 483     if ( !RewindReaders() ) {
 484         cerr << "BamMultiReader ERROR: could not rewind file(s) successfully";
 485         return false;
 486     }
 487
 488     // return status of cache update
 489     return UpdateAlignmentCache();
 490 }
 491
 492 // returns BAM file pointers to beginning of alignment data
 493 bool BamMultiReaderPrivate::RewindReaders(void) {
 494
 495     bool result = true;
 496
 497     // iterate over readers
 498     vector<MergeItem>::iterator readerIter = m_readers.begin();
 499     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 500     for ( ; readerIter != readerEnd; ++readerIter ) {
 501         MergeItem& item = (*readerIter);
 502         BamReader* reader = item.Reader;
 503         if ( reader == 0 ) continue;
 504
 505         // attempt rewind on BamReader
 506         result &= reader->Rewind();
 507     }
 508
 509     return result;
 510 }
 511
 512 void BamMultiReaderPrivate::SaveNextAlignment(BamReader* reader, BamAlignment* alignment) {
 513
 514     // if can read alignment from reader, store in cache
 515     // N.B. - lazy building of alignment's char data,
 516     // only populated on demand by sorting merger or client call to GetNextAlignment()
 517     if ( reader->GetNextAlignmentCore(*alignment) )
 518         m_alignmentCache->Add(MergeItem(reader, alignment));
 519 }
 520
 521 // sets the index caching mode on the readers
 522 void BamMultiReaderPrivate::SetIndexCacheMode(const BamIndex::IndexCacheMode mode) {
 523
 524     // iterate over readers
 525     vector<MergeItem>::iterator readerIter = m_readers.begin();
 526     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 527     for ( ; readerIter != readerEnd; ++readerIter ) {
 528         MergeItem& item = (*readerIter);
 529         BamReader* reader = item.Reader;
 530         if ( reader == 0 ) continue;
 531
 532         // set reader's index cache mode
 533         reader->SetIndexCacheMode(mode);
 534     }
 535 }
 536
 537 bool BamMultiReaderPrivate::SetRegion(const BamRegion& region) {
 538
 539     // NB: While it may make sense to track readers in which we can
 540     // successfully SetRegion, In practice a failure of SetRegion means "no
 541     // alignments here."  It makes sense to simply accept the failure,
 542     // UpdateAlignments(), and continue.
 543
 544     // iterate over alignments
 545     vector<MergeItem>::iterator readerIter = m_readers.begin();
 546     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 547     for ( ; readerIter != readerEnd; ++readerIter ) {
 548         MergeItem& item = (*readerIter);
 549         BamReader* reader = item.Reader;
 550         if ( reader == 0 ) continue;
 551
 552         // attempt to set BamReader's region of interest
 553         if ( !reader->SetRegion(region) ) {
 554             cerr << "BamMultiReader WARNING: could not jump " << reader->GetFilename() << " to "
 555                  << region.LeftRefID  << ":" << region.LeftPosition   << ".."
 556                  << region.RightRefID << ":" << region.RightPosition  << endl;
 557         }
 558     }
 559
 560     // return status of cache update
 561     return UpdateAlignmentCache();
 562 }
 563
 564 // updates our alignment cache
 565 bool BamMultiReaderPrivate::UpdateAlignmentCache(void) {
 566
 567     // create alignment cache if not created yet
 568     if ( m_alignmentCache == 0 ) {
 569         m_alignmentCache = CreateAlignmentCache();
 570         if ( m_alignmentCache == 0 ) {
 571             // set error string
 572             return false;
 573         }
 574     }
 575
 576     // clear any prior cache data
 577     m_alignmentCache->Clear();
 578
 579     // iterate over readers
 580     vector<MergeItem>::iterator readerIter = m_readers.begin();
 581     vector<MergeItem>::iterator readerEnd  = m_readers.end();
 582     for ( ; readerIter != readerEnd; ++readerIter ) {
 583         MergeItem& item = (*readerIter);
 584         BamReader* reader = item.Reader;
 585         BamAlignment* alignment = item.Alignment;
 586         if ( reader == 0 || alignment == 0 ) continue;
 587
 588         // save next alignment from each reader in cache
 589         SaveNextAlignment(reader, alignment);
 590     }
 591
 592     // if we get here, ok
 593     return true;
 594 }
 595
 596 // ValidateReaders checks that all the readers point to BAM files representing
 597 // alignments against the same set of reference sequences, and that the
 598 // sequences are identically ordered.  If these checks fail the operation of
 599 // the multireader is undefined, so we force program exit.
 600 bool BamMultiReaderPrivate::ValidateReaders(void) const {
 601
 602     // skip if no readers opened
 603     if ( m_readers.empty() )
 604         return true;
 605
 606     // retrieve first reader
 607     const MergeItem& firstItem = m_readers.front();
 608     const BamReader* firstReader = firstItem.Reader;
 609     if ( firstReader == 0 ) return false;
 610
 611     // retrieve first reader's header data
 612     const SamHeader& firstReaderHeader = firstReader->GetHeader();
 613     const string& firstReaderSortOrder = firstReaderHeader.SortOrder;
 614
 615     // retrieve first reader's reference data
 616     const RefVector& firstReaderRefData = firstReader->GetReferenceData();
 617     const int firstReaderRefCount = firstReader->GetReferenceCount();
 618     const int firstReaderRefSize = firstReaderRefData.size();
 619
 620     // iterate over all readers
 621     vector<MergeItem>::const_iterator readerIter = m_readers.begin();
 622     vector<MergeItem>::const_iterator readerEnd  = m_readers.end();
 623     for ( ; readerIter != readerEnd; ++readerIter ) {
 624         const MergeItem& item = (*readerIter);
 625         BamReader* reader = item.Reader;
 626         if ( reader == 0 ) continue;
 627
 628         // get current reader's header data
 629         const SamHeader& currentReaderHeader = reader->GetHeader();
 630         const string& currentReaderSortOrder = currentReaderHeader.SortOrder;
 631
 632         // check compatible sort order
 633         if ( currentReaderSortOrder != firstReaderSortOrder ) {
 634             // error string
 635             cerr << "BamMultiReader ERROR: mismatched sort order in " << reader->GetFilename()
 636                  << ", expected "  << firstReaderSortOrder
 637                  << ", but found " << currentReaderSortOrder << endl;
 638             return false;
 639         }
 640
 641         // get current reader's reference data
 642         const RefVector currentReaderRefData = reader->GetReferenceData();
 643         const int currentReaderRefCount = reader->GetReferenceCount();
 644         const int currentReaderRefSize  = currentReaderRefData.size();
 645
 646         // init reference data iterators
 647         RefVector::const_iterator firstRefIter   = firstReaderRefData.begin();
 648         RefVector::const_iterator firstRefEnd    = firstReaderRefData.end();
 649         RefVector::const_iterator currentRefIter = currentReaderRefData.begin();
 650
 651         // compare reference counts from BamReader ( & container size, in case of BR error)
 652         if ( (currentReaderRefCount != firstReaderRefCount) ||
 653              (firstReaderRefSize    != currentReaderRefSize) )
 654         {
 655             cerr << "BamMultiReader ERROR: mismatched number of references in " << reader->GetFilename()
 656                  << " expected " << firstReaderRefCount
 657                  << " reference sequences but only found " << currentReaderRefCount << endl;
 658             return false;
 659         }
 660
 661         // this will be ok; we just checked above that we have identically-sized sets of references
 662         // here we simply check if they are all, in fact, equal in content
 663         while ( firstRefIter != firstRefEnd ) {
 664             const RefData& firstRef   = (*firstRefIter);
 665             const RefData& currentRef = (*currentRefIter);
 666
 667             // compare reference name & length
 668             if ( (firstRef.RefName   != currentRef.RefName) ||
 669                  (firstRef.RefLength != currentRef.RefLength) )
 670             {
 671                 cerr << "BamMultiReader ERROR: mismatched references found in " << reader->GetFilename()
 672                      << " expected: " << endl;
 673
 674                 // print first reader's reference data
 675                 RefVector::const_iterator refIter = firstReaderRefData.begin();
 676                 RefVector::const_iterator refEnd  = firstReaderRefData.end();
 677                 for ( ; refIter != refEnd; ++refIter ) {
 678                     const RefData& entry = (*refIter);
 679                     cerr << entry.RefName << " " << entry.RefLength << endl;
 680                 }
 681
 682                 cerr << "but found: " << endl;
 683
 684                 // print current reader's reference data
 685                 refIter = currentReaderRefData.begin();
 686                 refEnd  = currentReaderRefData.end();
 687                 for ( ; refIter != refEnd; ++refIter ) {
 688                     const RefData& entry = (*refIter);
 689                     cerr << entry.RefName << " " << entry.RefLength << endl;
 690                 }
 691
 692                 return false;
 693             }
 694
 695             // update iterators
 696             ++firstRefIter;
 697             ++currentRefIter;
 698         }
 699     }
 700
 701     // if we get here, everything checks out
 702     return true;
 703 }