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