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: 23 December 2010 (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_isSortedByPosition(true)
  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 the BAM files
  44 void BamMultiReaderPrivate::Close(void) {
  45
  46     // clear out alignment cache
  47     m_alignments->Clear();
  48
  49     // iterate over readers
  50     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
  51     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
  52     for ( ; readerIter != readerEnd; ++readerIter ) {
  53
  54         // close reader
  55         BamReader*    reader    = (*readerIter).first;
  56         BamAlignment* alignment = (*readerIter).second;
  57         if ( reader ) reader->Close();
  58
  59         // delete pointers
  60         delete reader;
  61         reader = 0;
  62         delete alignment;
  63         alignment = 0;
  64     }
  65
  66     // clear out readers
  67     m_readers.clear();
  68
  69     // reset flags
  70     m_isCoreMode = false;
  71     m_isSortedByPosition = true;
  72 }
  73
  74 // saves index data to BAM index files (".bai"/".bti") where necessary, returns success/fail
  75 bool BamMultiReaderPrivate::CreateIndexes(bool useStandardIndex) {
  76
  77     bool result = true;
  78     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
  79     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
  80     for ( ; readerIter != readerEnd; ++readerIter ) {
  81         BamReader* reader = (*readerIter).first;
  82         result &= reader->CreateIndex(useStandardIndex);
  83     }
  84     return result;
  85 }
  86
  87 const string BamMultiReaderPrivate::ExtractReadGroup(const string& headerLine) const {
  88
  89     string readGroup("");
  90     stringstream headerLineSs(headerLine);
  91     string part;
  92
  93     // parse @RG header line, looking for the ID: tag
  94     while( getline(headerLineSs, part, '\t') ) {
  95         stringstream partSs(part);
  96         string subtag;
  97         getline(partSs, subtag, ':');
  98         if ( subtag == "ID" ) {
  99             getline(partSs, readGroup, ':');
 100             break;
 101         }
 102     }
 103     return readGroup;
 104 }
 105
 106 // makes a virtual, unified header for all the bam files in the multireader
 107 const string BamMultiReaderPrivate::GetHeaderText(void) const {
 108
 109     // just spit single header out if only have one reader open
 110     if ( m_readers.size() == 1 ) {
 111
 112
 113         vector<ReaderAlignment>::const_iterator readerBegin = m_readers.begin();
 114         const ReaderAlignment& entry = (*readerBegin);
 115         const BamReader* reader = entry.first;
 116         if ( reader == 0 ) return "";
 117         return reader->GetHeaderText();
 118     }
 119
 120     string mergedHeader = "";
 121     map<string, bool> readGroups;
 122
 123     // foreach extraction entry (each BAM file)
 124     vector<ReaderAlignment>::const_iterator readerBegin = m_readers.begin();
 125     vector<ReaderAlignment>::const_iterator readerIter  = readerBegin;
 126     vector<ReaderAlignment>::const_iterator readerEnd   = m_readers.end();
 127     for ( ; readerIter != readerEnd; ++readerIter ) {
 128
 129         // get header from reader
 130         const BamReader* reader = (*readerIter).first;
 131         if ( reader == 0 ) continue;
 132         string headerText = reader->GetHeaderText();
 133         if ( headerText.empty() ) continue;
 134
 135         // store header text in lines
 136         map<string, bool> currentFileReadGroups;
 137         const vector<string> lines = SplitHeaderText(headerText);
 138
 139         // iterate over header lines
 140         vector<string>::const_iterator linesIter = lines.begin();
 141         vector<string>::const_iterator linesEnd  = lines.end();
 142         for ( ; linesIter != linesEnd; ++linesIter ) {
 143
 144             // get next line from header, skip if empty
 145             const string headerLine = (*linesIter);
 146             if ( headerLine.empty() ) { continue; }
 147
 148             // if first file, save HD & SQ entries
 149             if ( readerIter == readerBegin ) {
 150                 if ( headerLine.find("@HD") == 0 || headerLine.find("@SQ") == 0) {
 151                     mergedHeader.append(headerLine.c_str());
 152                     mergedHeader.append(1, '\n');
 153                 }
 154             }
 155
 156             // (for all files) append RG entries if they are unique
 157             if ( headerLine.find("@RG") == 0 ) {
 158
 159                 // extract read group name from line
 160                 const string readGroup = ExtractReadGroup(headerLine);
 161
 162                 // make sure not to duplicate @RG entries
 163                 if ( readGroups.find(readGroup) == readGroups.end() ) {
 164                     mergedHeader.append(headerLine.c_str() );
 165                     mergedHeader.append(1, '\n');
 166                     readGroups[readGroup] = true;
 167                     currentFileReadGroups[readGroup] = true;
 168                 } else {
 169                     // warn iff we are reading one file and discover duplicated @RG tags in the header
 170                     // otherwise, we emit no warning, as we might be merging multiple BAM files with identical @RG tags
 171                     if ( currentFileReadGroups.find(readGroup) != currentFileReadGroups.end() ) {
 172                         cerr << "WARNING: duplicate @RG tag " << readGroup
 173                             << " entry in header of " << reader->GetFilename() << endl;
 174                     }
 175                 }
 176             }
 177         }
 178     }
 179
 180     // return merged header text
 181     return mergedHeader;
 182 }
 183
 184 // get next alignment among all files
 185 bool BamMultiReaderPrivate::GetNextAlignment(BamAlignment& al) {
 186     return LoadNextAlignment(al, false);
 187 }
 188
 189 // get next alignment among all files without parsing character data from alignments
 190 bool BamMultiReaderPrivate::GetNextAlignmentCore(BamAlignment& al) {
 191     return LoadNextAlignment(al, true);
 192 }
 193
 194 // ---------------------------------------------------------------------------------------
 195 //
 196 // NB: The following GetReferenceX() functions assume that we have identical
 197 // references for all BAM files.  We enforce this by invoking the above
 198 // validation function (ValidateReaders) to verify that our reference data
 199 // is the same across all files on Open, so we will not encounter a situation
 200 // in which there is a mismatch and we are still live.
 201 //
 202 // ---------------------------------------------------------------------------------------
 203
 204 // returns the number of reference sequences
 205 const int BamMultiReaderPrivate::GetReferenceCount(void) const {
 206     const ReaderAlignment& firstReader = m_readers.front();
 207     const BamReader* reader = firstReader.first;
 208     if ( reader ) return reader->GetReferenceCount();
 209     else return 0;
 210 }
 211
 212 // returns vector of reference objects
 213 const RefVector BamMultiReaderPrivate::GetReferenceData(void) const {
 214     const ReaderAlignment& firstReader = m_readers.front();
 215     const BamReader* reader = firstReader.first;
 216     if ( reader ) return reader->GetReferenceData();
 217     else return RefVector();
 218 }
 219
 220 // returns refID from reference name
 221 const int BamMultiReaderPrivate::GetReferenceID(const string& refName) const {
 222     const ReaderAlignment& firstReader = m_readers.front();
 223     const BamReader* reader = firstReader.first;
 224     if ( reader ) return reader->GetReferenceID(refName);
 225     else return -1; // ERROR case - how to report
 226 }
 227
 228 // ---------------------------------------------------------------------------------------
 229
 230 // checks if any readers still have alignments
 231 bool BamMultiReaderPrivate::HasOpenReaders(void) {
 232     return ( m_alignments->Size() > 0 );
 233 }
 234
 235 // returns whether underlying BAM readers ALL have an index loaded
 236 // this is useful to indicate whether Jump() or SetRegion() are possible
 237 bool BamMultiReaderPrivate::IsIndexLoaded(void) const {
 238     bool ok = true;
 239     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 240     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 241     for ( ; readerIter != readerEnd; ++readerIter ) {
 242         const BamReader* reader = (*readerIter).first;
 243         if ( reader ) ok &= reader->IsIndexLoaded();
 244     }
 245     return ok;
 246 }
 247
 248 // jumps to specified region(refID, leftBound) in BAM files, returns success/fail
 249 bool BamMultiReaderPrivate::Jump(int refID, int position) {
 250
 251     bool ok = true;
 252     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 253     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 254     for ( ; readerIter != readerEnd; ++readerIter ) {
 255         BamReader* reader = (*readerIter).first;
 256         if ( reader == 0 ) continue;
 257
 258         ok &= reader->Jump(refID, position);
 259         if ( !ok ) {
 260             cerr << "ERROR: could not jump " << reader->GetFilename()
 261                  << " to " << refID << ":" << position << endl;
 262             exit(1);
 263         }
 264     }
 265
 266     if (ok) UpdateAlignments();
 267     return ok;
 268 }
 269
 270 bool BamMultiReaderPrivate::LoadNextAlignment(BamAlignment& al, bool coreMode) {
 271
 272     // bail out if no more data to process
 273     if ( !HasOpenReaders() ) return false;
 274
 275     // "pop" next alignment and reader
 276     ReaderAlignment nextReaderAlignment = m_alignments->TakeFirst();
 277     BamReader*    reader    = nextReaderAlignment.first;
 278     BamAlignment* alignment = nextReaderAlignment.second;
 279
 280     // save it by copy to our argument
 281     al = BamAlignment(*alignment);
 282
 283     // peek to next alignment & store in cache
 284     m_isCoreMode = coreMode;
 285     SaveNextAlignment(reader,alignment);
 286
 287     // return success
 288     return true;
 289 }
 290
 291 // opens BAM files
 292 bool BamMultiReaderPrivate::Open(const vector<string>& filenames,
 293                                  bool openIndexes,
 294                                  bool coreMode,
 295                                  bool preferStandardIndex)
 296 {
 297     // store core mode flag
 298     m_isCoreMode = coreMode;
 299
 300     // first clear out any prior alignment cache prior data
 301     if ( m_alignments ) {
 302         m_alignments->Clear();
 303         delete m_alignments;
 304         m_alignments = 0;
 305     }
 306
 307     // create alignment cache based on sorting mode
 308     if ( m_isSortedByPosition )
 309         m_alignments = new PositionMultiMerger;
 310     else
 311         m_alignments = new ReadNameMultiMerger;
 312
 313     // iterate over filenames
 314     vector<string>::const_iterator filenameIter = filenames.begin();
 315     vector<string>::const_iterator filenameEnd  = filenames.end();
 316     for ( ; filenameIter != filenameEnd; ++filenameIter ) {
 317         const string filename = (*filenameIter);
 318
 319         bool openedOk = true;
 320         BamReader* reader = new BamReader;
 321         openedOk = reader->Open(filename, "", openIndexes, preferStandardIndex);
 322
 323         // if file opened ok
 324         if ( openedOk ) {
 325
 326             // try to read first alignment
 327             bool fileOk = true;
 328             BamAlignment* alignment = new BamAlignment;
 329             fileOk &= ( coreMode ? reader->GetNextAlignmentCore(*alignment)
 330                                  : reader->GetNextAlignment(*alignment) );
 331
 332             if ( fileOk ) {
 333
 334                 m_readers.push_back( make_pair(reader, alignment) );
 335                 m_alignments->Add( make_pair(reader, alignment) );
 336
 337             } else {
 338                 cerr << "WARNING: could not read first alignment in "
 339                      << filename << ", ignoring file" << endl;
 340
 341                 // if only file available & could not be read, return failure
 342                 if ( filenames.size() == 1 )
 343                     return false;
 344             }
 345         }
 346
 347         // TODO; any further error handling when openedOK is false ??
 348         else return false;
 349     }
 350
 351     // files opened ok, at least one alignment could be read,
 352     // now need to check that all files use same reference data
 353     ValidateReaders();
 354     return true;
 355 }
 356
 357 // print associated filenames to stdout
 358 void BamMultiReaderPrivate::PrintFilenames(void) const {
 359
 360     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 361     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 362     for ( ; readerIter != readerEnd; ++readerIter ) {
 363         const BamReader* reader = (*readerIter).first;
 364         if ( reader == 0 ) continue;
 365         cout << reader->GetFilename() << endl;
 366     }
 367 }
 368
 369 // returns BAM file pointers to beginning of alignment data
 370 bool BamMultiReaderPrivate::Rewind(void) {
 371
 372     bool result = true;
 373     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 374     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 375     for ( ; readerIter != readerEnd; ++readerIter ) {
 376         BamReader* reader = (*readerIter).first;
 377         if ( reader == 0 ) continue;
 378         result &= reader->Rewind();
 379     }
 380     return result;
 381 }
 382
 383 void BamMultiReaderPrivate::SaveNextAlignment(BamReader* reader, BamAlignment* alignment) {
 384
 385     // must be in core mode && sorting by position to call GNACore()
 386     if ( m_isCoreMode && m_isSortedByPosition ) {
 387         if ( reader->GetNextAlignmentCore(*alignment) )
 388             m_alignments->Add( make_pair(reader, alignment) );
 389     }
 390
 391     // not in core mode and/or sorting by readname, must call GNA()
 392     else {
 393         if ( reader->GetNextAlignment(*alignment) )
 394             m_alignments->Add( make_pair(reader, alignment) );
 395     }
 396 }
 397
 398 // sets the index caching mode on the readers
 399 void BamMultiReaderPrivate::SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode) {
 400
 401     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 402     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 403     for ( ; readerIter != readerEnd; ++readerIter ) {
 404         BamReader* reader = (*readerIter).first;
 405         if ( reader == 0 ) continue;
 406         reader->SetIndexCacheMode(mode);
 407     }
 408 }
 409
 410 bool BamMultiReaderPrivate::SetRegion(const BamRegion& region) {
 411
 412     // NB: While it may make sense to track readers in which we can
 413     // successfully SetRegion, In practice a failure of SetRegion means "no
 414     // alignments here."  It makes sense to simply accept the failure,
 415     // UpdateAlignments(), and continue.
 416
 417     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 418     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 419     for ( ; readerIter != readerEnd; ++readerIter ) {
 420         BamReader* reader = (*readerIter).first;
 421         if ( reader == 0 ) continue;
 422         if ( !reader->SetRegion(region) ) {
 423             cerr << "ERROR: could not jump " << reader->GetFilename() << " to "
 424                  << region.LeftRefID  << ":" << region.LeftPosition   << ".."
 425                  << region.RightRefID << ":" << region.RightPosition  << endl;
 426         }
 427     }
 428
 429     UpdateAlignments();
 430     return true;
 431 }
 432
 433 void BamMultiReaderPrivate::SetSortOrder(const BamMultiReader::SortOrder& order) {
 434
 435     const BamMultiReader::SortOrder oldSortOrder = ( m_isSortedByPosition ? BamMultiReader::SortedByPosition
 436                                                                           : BamMultiReader::SortedByReadName ) ;
 437     // skip if no change needed
 438     if ( oldSortOrder == order ) return;
 439
 440     // create new alignment cache
 441     IBamMultiMerger* newAlignmentCache(0);
 442     if ( order == BamMultiReader::SortedByPosition ) {
 443         newAlignmentCache = new PositionMultiMerger;
 444         m_isSortedByPosition = true;
 445     }
 446     else {
 447         newAlignmentCache = new ReadNameMultiMerger;
 448         m_isSortedByPosition = false;
 449     }
 450
 451     // copy old cache contents to new cache
 452     while ( m_alignments->Size() > 0 ) {
 453         ReaderAlignment value = m_alignments->TakeFirst();
 454         newAlignmentCache->Add(value);
 455     }
 456
 457     // remove old cache structure & point to new cache
 458     delete m_alignments;
 459     m_alignments = newAlignmentCache;
 460 }
 461
 462 // updates our alignment cache
 463 void BamMultiReaderPrivate::UpdateAlignments(void) {
 464
 465     // clear the cache
 466     m_alignments->Clear();
 467
 468     // iterate over readers
 469     vector<ReaderAlignment>::iterator readerIter = m_readers.begin();
 470     vector<ReaderAlignment>::iterator readerEnd  = m_readers.end();
 471     for ( ; readerIter != readerEnd; ++readerIter ) {
 472         BamReader* reader = (*readerIter).first;
 473         BamAlignment* alignment = (*readerIter).second;
 474         if ( reader == 0 ) continue;
 475         SaveNextAlignment(reader, alignment);
 476     }
 477 }
 478
 479 // splits the entire header into a list of strings
 480 const vector<string> BamMultiReaderPrivate::SplitHeaderText(const string& headerText) const {
 481     stringstream header(headerText);
 482     vector<string> lines;
 483     string item;
 484     while ( getline(header, item) )
 485         lines.push_back(item);
 486     return lines;
 487 }
 488
 489 // ValidateReaders checks that all the readers point to BAM files representing
 490 // alignments against the same set of reference sequences, and that the
 491 // sequences are identically ordered.  If these checks fail the operation of
 492 // the multireader is undefined, so we force program exit.
 493 void BamMultiReaderPrivate::ValidateReaders(void) const {
 494
 495     // retrieve first reader data
 496     const BamReader* firstReader = m_readers.front().first;
 497     if ( firstReader == 0 ) return; // signal error?
 498     const RefVector firstReaderRefData = firstReader->GetReferenceData();
 499     const int firstReaderRefCount = firstReader->GetReferenceCount();
 500     const int firstReaderRefSize = firstReaderRefData.size();
 501
 502     // iterate over all readers
 503     vector<ReaderAlignment>::const_iterator readerIter = m_readers.begin();
 504     vector<ReaderAlignment>::const_iterator readerEnd  = m_readers.end();
 505     for ( ; readerIter != readerEnd; ++readerIter ) {
 506
 507         // get current reader data
 508         BamReader* reader = (*readerIter).first;
 509         if ( reader == 0 ) continue; // error?
 510         const RefVector currentReaderRefData = reader->GetReferenceData();
 511         const int currentReaderRefCount = reader->GetReferenceCount();
 512         const int currentReaderRefSize  = currentReaderRefData.size();
 513
 514         // init container iterators
 515         RefVector::const_iterator firstRefIter   = firstReaderRefData.begin();
 516         RefVector::const_iterator firstRefEnd    = firstReaderRefData.end();
 517         RefVector::const_iterator currentRefIter = currentReaderRefData.begin();
 518
 519         // compare reference counts from BamReader ( & container size, in case of BR error)
 520         if ( (currentReaderRefCount != firstReaderRefCount) ||
 521              (firstReaderRefSize    != currentReaderRefSize) )
 522         {
 523             cerr << "ERROR: mismatched number of references in " << reader->GetFilename()
 524                  << " expected " << firstReaderRefCount
 525                  << " reference sequences but only found " << currentReaderRefCount << endl;
 526             exit(1);
 527         }
 528
 529         // this will be ok; we just checked above that we have identically-sized sets of references
 530         // here we simply check if they are all, in fact, equal in content
 531         while ( firstRefIter != firstRefEnd ) {
 532
 533             const RefData& firstRef   = (*firstRefIter);
 534             const RefData& currentRef = (*currentRefIter);
 535
 536             // compare reference name & length
 537             if ( (firstRef.RefName   != currentRef.RefName) ||
 538                  (firstRef.RefLength != currentRef.RefLength) )
 539             {
 540                 cerr << "ERROR: mismatched references found in " << reader->GetFilename()
 541                      << " expected: " << endl;
 542                 RefVector::const_iterator refIter = firstReaderRefData.begin();
 543                 RefVector::const_iterator refEnd  = firstReaderRefData.end();
 544                 for ( ; refIter != refEnd; ++refIter ) {
 545                     const RefData& entry = (*refIter);
 546                     cerr << entry.RefName << " " << entry.RefLength << endl;
 547                 }
 548
 549                 cerr << "but found: " << endl;
 550                 refIter = currentReaderRefData.begin();
 551                 refEnd  = currentReaderRefData.end();
 552                 for ( ; refIter != refEnd; ++refIter ) {
 553                     const RefData& entry = (*refIter);
 554                     cerr << entry.RefName << " " << entry.RefLength << endl;
 555                 }
 556
 557                 exit(1);
 558             }
 559
 560             // update iterators
 561             ++firstRefIter;
 562             ++currentRefIter;
 563         }
 564     }
 565 }