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 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 9 September 2011 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Functionality for simultaneously reading multiple BAM files
   9 // *************************************************************************
  10
  11 #include <api/BamAlignment.h>
  12 #include <api/BamMultiReader.h>
  13 #include <api/internal/BamMultiMerger_p.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_alignments(0)
  28     , m_sortOrder(BamMultiReader::SortedByPosition)
  29 { }
  30
  31 // dtor
  32 BamMultiReaderPrivate::~BamMultiReaderPrivate(void) {
  33
  34     // close all open BAM readers
  35     Close();
  36
  37     // clean up alignment cache
  38     delete m_alignments;
  39     m_alignments = 0;
  40 }
  41
  42 // close all BAM files
  43 void BamMultiReaderPrivate::Close(void) {
  44     CloseFiles( Filenames() );
  45 }
  46
  47 // close requested BAM file
  48 void BamMultiReaderPrivate::CloseFile(const string& filename) {
  49     vector<string> filenames(1, filename);
  50     CloseFiles(filenames);
  51 }
  52
  53 // close requested BAM files
  54 void BamMultiReaderPrivate::CloseFiles(const vector<string>& filenames) {
  55
  56     // iterate over filenames
  57     vector<string>::const_iterator filesIter = filenames.begin();
  58     vector<string>::const_iterator filesEnd  = filenames.end();
  59     for ( ; filesIter != filesEnd; ++filesIter ) {
  60         const string& filename = (*filesIter);
  61         if ( filename.empty() ) continue;
  62
  63         // iterate over readers
  64         vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
  65         vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
  66         for ( ; readerIter != readerEnd; ++readerIter ) {
  67             BamReader* reader = (*readerIter).first;
  68             if ( reader == 0 ) continue;
  69
  70             // if reader matches requested filename
  71             if ( reader->GetFilename() == filename ) {
  72
  73                 // remove reader/alignment from alignment cache
  74                 m_alignments->Remove(reader);
  75
  76                 // close & delete reader
  77                 reader->Close();
  78                 delete reader;
  79                 reader = 0;
  80
  81                 // delete reader's alignment entry
  82                 BamAlignment* alignment = (*readerIter).second;
  83                 delete alignment;
  84                 alignment = 0;
  85
  86                 // remove reader from container
  87                 m_readers.erase(readerIter);
  88
  89                 // on match, just go on to next filename
  90                 // (no need to keep looking and iterator is invalid now anyway)
  91                 break;
  92             }
  93         }
  94     }
  95
  96     // make sure alignment cache is cleared if all readers are now closed
  97     if ( m_readers.empty() && m_alignments != 0 )
  98         m_alignments->Clear();
  99 }
 100
 101 // creates index files for BAM files that don't have them
 102 bool BamMultiReaderPrivate::CreateIndexes(const BamIndex::IndexType& type) {
 103
 104     bool result = true;
 105
 106     // iterate over readers
 107     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 108     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 109     for ( ; readerIter != readerEnd; ++readerIter ) {
 110         BamReader* reader = (*readerIter).first;
 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 IBamMultiMerger* BamMultiReaderPrivate::CreateMergerForCurrentSortOrder(void) const {
 122     switch ( m_sortOrder ) {
 123         case ( BamMultiReader::SortedByPosition ) : return new PositionMultiMerger;
 124         case ( BamMultiReader::SortedByReadName ) : return new ReadNameMultiMerger;
 125         case ( BamMultiReader::Unsorted )         : return new UnsortedMultiMerger;
 126         default :
 127             cerr << "BamMultiReader ERROR: requested sort order is unknown" << endl;
 128             return 0;
 129     }
 130 }
 131
 132 const string BamMultiReaderPrivate::ExtractReadGroup(const string& headerLine) const {
 133
 134     string readGroup("");
 135     stringstream headerLineSs(headerLine);
 136     string part;
 137
 138     // parse @RG header line, looking for the ID: tag
 139     while( getline(headerLineSs, part, '\t') ) {
 140         stringstream partSs(part);
 141         string subtag;
 142         getline(partSs, subtag, ':');
 143         if ( subtag == "ID" ) {
 144             getline(partSs, readGroup, ':');
 145             break;
 146         }
 147     }
 148     return readGroup;
 149 }
 150
 151 const vector<string> BamMultiReaderPrivate::Filenames(void) const {
 152
 153     // init filename container
 154     vector<string> filenames;
 155     filenames.reserve( m_readers.size() );
 156
 157     // iterate over readers
 158     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 159     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 160     for ( ; readerIter != readerEnd; ++readerIter ) {
 161         const BamReader* reader = (*readerIter).first;
 162         if ( reader == 0 ) continue;
 163
 164         // store filename if not empty
 165         const string filename = reader->GetFilename();
 166         if ( !filename.empty() )
 167             filenames.push_back( reader->GetFilename() );
 168     }
 169
 170     // return result
 171     return filenames;
 172 }
 173
 174 SamHeader BamMultiReaderPrivate::GetHeader(void) const {
 175     string text = GetHeaderText();
 176     return SamHeader(text);
 177 }
 178
 179 // makes a virtual, unified header for all the bam files in the multireader
 180 string BamMultiReaderPrivate::GetHeaderText(void) const {
 181
 182     // TODO: merge SamHeader objects instead of parsing string data (again)
 183
 184     // if only one reader is open
 185     if ( m_readers.size() == 1 ) {
 186
 187         // just return reader's header text
 188         const ReaderAlignment& ra = m_readers.front();
 189         const BamReader* reader = ra.first;
 190         if ( reader ) return reader->GetHeaderText();
 191
 192         // invalid reader
 193         return string();
 194     }
 195
 196     string mergedHeader("");
 197     map<string, bool> readGroups;
 198
 199     // foreach extraction entry (each BAM file)
 200     vector<ReaderAlignment>::const_iterator readerBegin = m_readers.begin();
 201     vector<ReaderAlignment>::const_iterator readerIter  = readerBegin;
 202     vector<ReaderAlignment>::const_iterator readerEnd   = m_readers.end();
 203     for ( ; readerIter != readerEnd; ++readerIter ) {
 204         const BamReader* reader = (*readerIter).first;
 205         if ( reader == 0 ) continue;
 206
 207         // get header from reader
 208         string headerText = reader->GetHeaderText();
 209         if ( headerText.empty() ) continue;
 210
 211         // store header text in lines
 212         map<string, bool> currentFileReadGroups;
 213         const vector<string> lines = SplitHeaderText(headerText);
 214
 215         // iterate over header lines
 216         vector<string>::const_iterator linesIter = lines.begin();
 217         vector<string>::const_iterator linesEnd  = lines.end();
 218         for ( ; linesIter != linesEnd; ++linesIter ) {
 219
 220             // get next line from header, skip if empty
 221             const string headerLine = (*linesIter);
 222             if ( headerLine.empty() ) continue;
 223
 224             // if first file, save HD & SQ entries
 225             // TODO: what if first file has empty header, should just check for empty 'mergedHeader' instead ?
 226             if ( readerIter == readerBegin ) {
 227                 if ( headerLine.find("@HD") == 0 || headerLine.find("@SQ") == 0) {
 228                     mergedHeader.append(headerLine.c_str());
 229                     mergedHeader.append(1, '\n');
 230                 }
 231             }
 232
 233             // (for all files) append RG entries if they are unique
 234             if ( headerLine.find("@RG") == 0 ) {
 235
 236                 // extract read group name from line
 237                 const string readGroup = ExtractReadGroup(headerLine);
 238
 239                 // make sure not to duplicate @RG entries
 240                 if ( readGroups.find(readGroup) == readGroups.end() ) {
 241                     mergedHeader.append(headerLine.c_str() );
 242                     mergedHeader.append(1, '\n');
 243                     readGroups[readGroup] = true;
 244                     currentFileReadGroups[readGroup] = true;
 245                 } else {
 246                     // warn iff we are reading one file and discover duplicated @RG tags in the header
 247                     // otherwise, we emit no warning, as we might be merging multiple BAM files with identical @RG tags
 248                     if ( currentFileReadGroups.find(readGroup) != currentFileReadGroups.end() ) {
 249                         cerr << "BamMultiReader WARNING: duplicate @RG tag " << readGroup
 250                              << " entry in header of " << reader->GetFilename() << endl;
 251                     }
 252                 }
 253             }
 254         }
 255     }
 256
 257     // return merged header text
 258     return mergedHeader;
 259 }
 260
 261 // get next alignment among all files
 262 bool BamMultiReaderPrivate::GetNextAlignment(BamAlignment& al) {
 263     return PopNextCachedAlignment(al, true);
 264 }
 265
 266 // get next alignment among all files without parsing character data from alignments
 267 bool BamMultiReaderPrivate::GetNextAlignmentCore(BamAlignment& al) {
 268     return PopNextCachedAlignment(al, false);
 269 }
 270
 271 // ---------------------------------------------------------------------------------------
 272 //
 273 // NB: The following GetReferenceX() functions assume that we have identical
 274 // references for all BAM files.  We enforce this by invoking the
 275 // ValidateReaders() method to verify that our reference data is the same
 276 // across all files on Open - so we will not encounter a situation in which
 277 // there is a mismatch and we are still live.
 278 //
 279 // ---------------------------------------------------------------------------------------
 280
 281 // returns the number of reference sequences
 282 int BamMultiReaderPrivate::GetReferenceCount(void) const {
 283
 284     // handle empty multireader
 285     if ( m_readers.empty() )
 286         return 0;
 287
 288     // return reference count from first reader
 289     const ReaderAlignment& ra = m_readers.front();
 290     const BamReader* reader = ra.first;
 291     if ( reader ) return reader->GetReferenceCount();
 292
 293     // invalid reader
 294     return 0;
 295 }
 296
 297 // returns vector of reference objects
 298 const RefVector BamMultiReaderPrivate::GetReferenceData(void) const {
 299
 300     // handle empty multireader
 301     if ( m_readers.empty() )
 302         return RefVector();
 303
 304     // return reference data from first BamReader
 305     const ReaderAlignment& ra = m_readers.front();
 306     const BamReader* reader = ra.first;
 307     if ( reader ) return reader->GetReferenceData();
 308
 309     // invalid reader
 310     return RefVector();
 311 }
 312
 313 // returns refID from reference name
 314 int BamMultiReaderPrivate::GetReferenceID(const string& refName) const {
 315
 316     // handle empty multireader
 317     if ( m_readers.empty() )
 318         return -1;
 319
 320     // return reference ID from first BamReader
 321     const ReaderAlignment& ra = m_readers.front();
 322     const BamReader* reader = ra.first;
 323     if ( reader ) return reader->GetReferenceID(refName);
 324
 325     // invalid reader
 326     return -1;
 327 }
 328 // ---------------------------------------------------------------------------------------
 329
 330 // checks if any readers still have alignments
 331 bool BamMultiReaderPrivate::HasAlignmentData(void) const {
 332     if ( m_alignments == 0 )
 333         return false;
 334     return !m_alignments->IsEmpty();
 335 }
 336
 337 // returns true if all readers have index data available
 338 // this is useful to indicate whether Jump() or SetRegion() are possible
 339 bool BamMultiReaderPrivate::HasIndexes(void) const {
 340
 341     // handle empty multireader
 342     if ( m_readers.empty() )
 343         return false;
 344
 345     bool result = true;
 346
 347     // iterate over readers
 348     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 349     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 350     for ( ; readerIter != readerEnd; ++readerIter ) {
 351         const BamReader* reader = (*readerIter).first;
 352         if ( reader  == 0 ) continue;
 353
 354         // see if current reader has index data
 355         result &= reader->HasIndex();
 356     }
 357
 358     return result;
 359 }
 360
 361 // returns true if multireader has open readers
 362 bool BamMultiReaderPrivate::HasOpenReaders(void) {
 363
 364     // iterate over readers
 365     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 366     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 367     for ( ; readerIter != readerEnd; ++readerIter ) {
 368         const BamReader* reader = (*readerIter).first;
 369         if ( reader == 0 ) continue;
 370
 371         // return true whenever an open reader is found
 372         if ( reader->IsOpen() ) return true;
 373     }
 374
 375     // no readers open
 376     return false;
 377 }
 378
 379 // performs random-access jump using (refID, position) as a left-bound
 380 bool BamMultiReaderPrivate::Jump(int refID, int position) {
 381
 382     // NB: While it may make sense to track readers in which we can
 383     // successfully Jump, in practice a failure of Jump means "no
 384     // alignments here."  It makes sense to simply accept the failure,
 385     // UpdateAlignments(), and continue.
 386
 387     // iterate over readers
 388     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 389     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 390     for ( ; readerIter != readerEnd; ++readerIter ) {
 391         BamReader* reader = (*readerIter).first;
 392         if ( reader == 0 ) continue;
 393
 394         // attempt jump() on each
 395         if ( !reader->Jump(refID, position) ) {
 396             cerr << "BamMultiReader ERROR: could not jump " << reader->GetFilename()
 397                  << " to " << refID << ":" << position << endl;
 398         }
 399     }
 400
 401     // update alignment cache & return success
 402     UpdateAlignmentCache();
 403     return true;
 404 }
 405
 406 bool BamMultiReaderPrivate::LoadNextAlignment(BamReader* reader, BamAlignment* alignment) {
 407     // lazy building of alignment's char data,
 408     // only populated on demand by sorting merger or client call to GetNextAlignment()
 409     return reader->GetNextAlignmentCore(*alignment);
 410 }
 411
 412 // locate (& load) index files for BAM readers that don't already have one loaded
 413 bool BamMultiReaderPrivate::LocateIndexes(const BamIndex::IndexType& preferredType) {
 414
 415     bool result = true;
 416
 417     // iterate over readers
 418     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 419     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 420     for ( ; readerIter != readerEnd; ++readerIter ) {
 421         BamReader* reader = (*readerIter).first;
 422         if ( reader == 0 ) continue;
 423
 424         // if reader has no index, try to locate one
 425         if ( !reader->HasIndex() )
 426             result &= reader->LocateIndex(preferredType);
 427     }
 428
 429     return result;
 430 }
 431
 432 // opens BAM files
 433 bool BamMultiReaderPrivate::Open(const vector<string>& filenames) {
 434
 435     // create alignment cache if neccessary
 436     if ( m_alignments == 0 ) {
 437         m_alignments = CreateMergerForCurrentSortOrder();
 438         if ( m_alignments == 0 ) return false;
 439     }
 440
 441     // iterate over filenames
 442     vector<string>::const_iterator filenameIter = filenames.begin();
 443     vector<string>::const_iterator filenameEnd  = filenames.end();
 444     for ( ; filenameIter != filenameEnd; ++filenameIter ) {
 445         const string& filename = (*filenameIter);
 446         if ( filename.empty() ) continue;
 447
 448         // attempt to open BamReader on filename
 449         BamReader* reader = OpenReader(filename);
 450         if ( reader == 0 ) continue;
 451
 452         // store reader with new alignment
 453         m_readers.push_back( make_pair(reader, new BamAlignment) );
 454     }
 455
 456     // validate & rewind any opened readers, also refreshes alignment cache
 457     if ( !m_readers.empty() ) {
 458         ValidateReaders();
 459         Rewind();
 460     }
 461
 462     // return success
 463     return true;
 464 }
 465
 466 bool BamMultiReaderPrivate::OpenFile(const std::string& filename) {
 467     vector<string> filenames(1, filename);
 468     return Open(filenames);
 469 }
 470
 471 bool BamMultiReaderPrivate::OpenIndexes(const vector<string>& indexFilenames) {
 472
 473     // TODO: This needs to be cleaner - should not assume same order.
 474     //       And either way, shouldn't start at first reader.  Should start at
 475     //       first reader without an index?
 476
 477     // make sure same number of index filenames as readers
 478     if ( m_readers.size() != indexFilenames.size() || !indexFilenames.empty() )
 479         return false;
 480
 481     // init result flag
 482     bool result = true;
 483
 484     // iterate over BamReaders
 485     vector<string>::const_iterator indexFilenameIter = indexFilenames.begin();
 486     vector<string>::const_iterator indexFilenameEnd  = indexFilenames.end();
 487     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 488     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 489     for ( ; readerIter != readerEnd; ++readerIter ) {
 490         BamReader* reader = (*readerIter).first;
 491
 492         // open index filename on reader
 493         if ( reader ) {
 494             const string& indexFilename = (*indexFilenameIter);
 495             result &= reader->OpenIndex(indexFilename);
 496         }
 497
 498         // increment filename iterator, skip if no more index files to open
 499         if ( ++indexFilenameIter == indexFilenameEnd )
 500             break;
 501     }
 502
 503     // TODO: validation ??
 504
 505     // return success/fail
 506     return result;
 507 }
 508
 509 BamReader* BamMultiReaderPrivate::OpenReader(const std::string& filename) {
 510
 511     // create new BamReader
 512     BamReader* reader = new BamReader;
 513
 514     // if reader opens OK
 515     if ( reader->Open(filename) ) {
 516
 517         // attempt to read first alignment (sanity check)
 518         // if ok, then return BamReader pointer
 519         BamAlignment al;
 520         if ( reader->GetNextAlignmentCore(al) )
 521             return reader;
 522
 523         // could not read alignment
 524         else {
 525             cerr << "BamMultiReader WARNING: Could not read first alignment from "
 526                  << filename << ", ignoring file" << endl;
 527         }
 528     }
 529
 530     // reader could not open
 531     else {
 532         cerr << "BamMultiReader WARNING: Could not open "
 533               << filename << ", ignoring file" << endl;
 534     }
 535
 536     // if we get here, there was a problem with this BAM file (opening or reading)
 537     // clean up memory allocation & return null pointer
 538     delete reader;
 539     return 0;
 540 }
 541
 542 // print associated filenames to stdout
 543 void BamMultiReaderPrivate::PrintFilenames(void) const {
 544     const vector<string>& filenames = Filenames();
 545     vector<string>::const_iterator filenameIter = filenames.begin();
 546     vector<string>::const_iterator filenameEnd  = filenames.end();
 547     for ( ; filenameIter != filenameEnd; ++filenameIter )
 548         cout << (*filenameIter) << endl;
 549 }
 550
 551 bool BamMultiReaderPrivate::PopNextCachedAlignment(BamAlignment& al, const bool needCharData) {
 552
 553     // bail out if no more data to process
 554     if ( !HasAlignmentData() )
 555         return false;
 556
 557     // "pop" next alignment and reader
 558     ReaderAlignment nextReaderAlignment = m_alignments->TakeFirst();
 559     BamReader* reader = nextReaderAlignment.first;
 560     BamAlignment* alignment = nextReaderAlignment.second;
 561
 562     // store cached alignment into destination parameter (by copy)
 563     al = *alignment;
 564
 565     // set char data if requested
 566     if ( needCharData ) {
 567         al.BuildCharData();
 568         al.Filename = reader->GetFilename();
 569     }
 570
 571     // peek to next alignment & store in cache
 572     SaveNextAlignment(reader, alignment);
 573
 574     // return success
 575     return true;
 576 }
 577
 578 // returns BAM file pointers to beginning of alignment data & resets alignment cache
 579 bool BamMultiReaderPrivate::Rewind(void) {
 580
 581     // clear out alignment cache
 582     m_alignments->Clear();
 583
 584     // attempt to rewind files
 585     if ( !RewindReaders() ) {
 586         cerr << "BamMultiReader ERROR: could not rewind file(s) successfully";
 587         return false;
 588     }
 589
 590     // reset cache & return success
 591     UpdateAlignmentCache();
 592     return true;
 593 }
 594
 595 // returns BAM file pointers to beginning of alignment data
 596 bool BamMultiReaderPrivate::RewindReaders(void) {
 597
 598     bool result = true;
 599
 600     // iterate over readers
 601     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 602     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 603     for ( ; readerIter != readerEnd; ++readerIter ) {
 604         BamReader* reader = (*readerIter).first;
 605         if ( reader == 0 ) continue;
 606
 607         // attempt rewind on BamReader
 608         result &= reader->Rewind();
 609     }
 610
 611     return result;
 612 }
 613
 614 void BamMultiReaderPrivate::SaveNextAlignment(BamReader* reader, BamAlignment* alignment) {
 615
 616     // if can read alignment from reader, store in cache
 617     if ( LoadNextAlignment(reader, alignment) )
 618         m_alignments->Add( make_pair(reader, alignment) );
 619 }
 620
 621 // sets the index caching mode on the readers
 622 void BamMultiReaderPrivate::SetIndexCacheMode(const BamIndex::IndexCacheMode mode) {
 623
 624     // iterate over readers
 625     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 626     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 627     for ( ; readerIter != readerEnd; ++readerIter ) {
 628         BamReader* reader = (*readerIter).first;
 629         if ( reader == 0 ) continue;
 630
 631         // set reader's index cache mode
 632         reader->SetIndexCacheMode(mode);
 633     }
 634 }
 635
 636 bool BamMultiReaderPrivate::SetRegion(const BamRegion& region) {
 637
 638     // NB: While it may make sense to track readers in which we can
 639     // successfully SetRegion, In practice a failure of SetRegion means "no
 640     // alignments here."  It makes sense to simply accept the failure,
 641     // UpdateAlignments(), and continue.
 642
 643     // iterate over alignments
 644     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 645     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 646     for ( ; readerIter != readerEnd; ++readerIter ) {
 647         BamReader* reader = (*readerIter).first;
 648         if ( reader == 0 ) continue;
 649
 650         // attempt to set BamReader's region of interest
 651         if ( !reader->SetRegion(region) ) {
 652             cerr << "BamMultiReader WARNING: could not jump " << reader->GetFilename() << " to "
 653                  << region.LeftRefID  << ":" << region.LeftPosition   << ".."
 654                  << region.RightRefID << ":" << region.RightPosition  << endl;
 655         }
 656     }
 657
 658     // update alignment cache & return success
 659     UpdateAlignmentCache();
 660     return true;
 661 }
 662
 663 void BamMultiReaderPrivate::SetSortOrder(const BamMultiReader::SortOrder& order) {
 664
 665     // skip if no change needed
 666     if ( m_sortOrder == order ) return;
 667
 668     // set new sort order
 669     m_sortOrder = order;
 670
 671     // create new alignment cache based on sort order
 672     IBamMultiMerger* newAlignmentCache = CreateMergerForCurrentSortOrder();
 673     if ( newAlignmentCache == 0 ) return; // print error?
 674
 675     // copy old cache contents to new cache
 676     while ( m_alignments->Size() > 0 ) {
 677         ReaderAlignment value = m_alignments->TakeFirst(); // retrieves & 'pops'
 678         newAlignmentCache->Add(value);
 679     }
 680
 681     // remove old cache structure & point to new cache
 682     delete m_alignments;
 683     m_alignments = newAlignmentCache;
 684 }
 685
 686 // splits the entire header into a list of strings
 687 const vector<string> BamMultiReaderPrivate::SplitHeaderText(const string& headerText) const {
 688
 689     stringstream header(headerText);
 690     string item;
 691
 692     vector<string> lines;
 693     while ( getline(header, item) )
 694         lines.push_back(item);
 695     return lines;
 696 }
 697
 698 // updates our alignment cache
 699 void BamMultiReaderPrivate::UpdateAlignmentCache(void) {
 700
 701     // skip if invalid alignment cache
 702     if ( m_alignments == 0 ) return;
 703
 704     // clear the cache
 705     m_alignments->Clear();
 706
 707     // iterate over readers
 708     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 709     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 710     for ( ; readerIter != readerEnd; ++readerIter ) {
 711         BamReader* reader = (*readerIter).first;
 712         BamAlignment* alignment = (*readerIter).second;
 713         if ( reader == 0 || alignment == 0 ) continue;
 714
 715         // save next alignment from each reader in cache
 716         SaveNextAlignment(reader, alignment);
 717     }
 718 }
 719
 720 // ValidateReaders checks that all the readers point to BAM files representing
 721 // alignments against the same set of reference sequences, and that the
 722 // sequences are identically ordered.  If these checks fail the operation of
 723 // the multireader is undefined, so we force program exit.
 724 void BamMultiReaderPrivate::ValidateReaders(void) const {
 725
 726     // retrieve first reader data
 727     const BamReader* firstReader = m_readers.front().first;
 728     if ( firstReader == 0 ) return;
 729     const RefVector firstReaderRefData = firstReader->GetReferenceData();
 730     const int firstReaderRefCount = firstReader->GetReferenceCount();
 731     const int firstReaderRefSize = firstReaderRefData.size();
 732
 733     // iterate over all readers
 734     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 735     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 736     for ( ; readerIter != readerEnd; ++readerIter ) {
 737
 738         // get current reader data
 739         BamReader* reader = (*readerIter).first;
 740         if ( reader == 0 ) continue;
 741         const RefVector currentReaderRefData = reader->GetReferenceData();
 742         const int currentReaderRefCount = reader->GetReferenceCount();
 743         const int currentReaderRefSize  = currentReaderRefData.size();
 744
 745         // init container iterators
 746         RefVector::const_iterator firstRefIter   = firstReaderRefData.begin();
 747         RefVector::const_iterator firstRefEnd    = firstReaderRefData.end();
 748         RefVector::const_iterator currentRefIter = currentReaderRefData.begin();
 749
 750         // compare reference counts from BamReader ( & container size, in case of BR error)
 751         if ( (currentReaderRefCount != firstReaderRefCount) ||
 752              (firstReaderRefSize    != currentReaderRefSize) )
 753         {
 754             cerr << "BamMultiReader ERROR: mismatched number of references in " << reader->GetFilename()
 755                  << " expected " << firstReaderRefCount
 756                  << " reference sequences but only found " << currentReaderRefCount << endl;
 757             exit(1);
 758         }
 759
 760         // this will be ok; we just checked above that we have identically-sized sets of references
 761         // here we simply check if they are all, in fact, equal in content
 762         while ( firstRefIter != firstRefEnd ) {
 763             const RefData& firstRef   = (*firstRefIter);
 764             const RefData& currentRef = (*currentRefIter);
 765
 766             // compare reference name & length
 767             if ( (firstRef.RefName   != currentRef.RefName) ||
 768                  (firstRef.RefLength != currentRef.RefLength) )
 769             {
 770                 cerr << "BamMultiReader ERROR: mismatched references found in " << reader->GetFilename()
 771                      << " expected: " << endl;
 772
 773                 // print first reader's reference data
 774                 RefVector::const_iterator refIter = firstReaderRefData.begin();
 775                 RefVector::const_iterator refEnd  = firstReaderRefData.end();
 776                 for ( ; refIter != refEnd; ++refIter ) {
 777                     const RefData& entry = (*refIter);
 778                     cerr << entry.RefName << " " << entry.RefLength << endl;
 779                 }
 780
 781                 cerr << "but found: " << endl;
 782
 783                 // print current reader's reference data
 784                 refIter = currentReaderRefData.begin();
 785                 refEnd  = currentReaderRefData.end();
 786                 for ( ; refIter != refEnd; ++refIter ) {
 787                     const RefData& entry = (*refIter);
 788                     cerr << entry.RefName << " " << entry.RefLength << endl;
 789                 }
 790
 791                 exit(1);
 792             }
 793
 794             // update iterators
 795             ++firstRefIter;
 796             ++currentRefIter;
 797         }
 798     }
 799 }