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