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