src/api/BamMultiReader.cpp

   1 // ***************************************************************************
   2 // BamMultiReader.cpp (c) 2010 Erik Garrison, Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 18 September 2010 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad
   9 // Institute.
  10 // ---------------------------------------------------------------------------
  11 // Functionality for simultaneously reading multiple BAM files.
  12 //
  13 // This functionality allows applications to work on very large sets of files
  14 // without requiring intermediate merge, sort, and index steps for each file
  15 // subset.  It also improves the performance of our merge system as it
  16 // precludes the need to sort merged files.
  17 // ***************************************************************************
  18
  19 #include <algorithm>
  20 #include <fstream>
  21 #include <iostream>
  22 #include <iterator>
  23 #include <sstream>
  24 #include <string>
  25 #include <vector>
  26 #include "BGZF.h"
  27 #include "BamMultiReader.h"
  28 using namespace BamTools;
  29 using namespace std;
  30
  31 // -----------------------------------------------------
  32 // BamMultiReader implementation
  33 // -----------------------------------------------------
  34
  35 // constructor
  36 BamMultiReader::BamMultiReader(void)
  37     : CurrentRefID(0)
  38     , CurrentLeft(0)
  39 { }
  40
  41 // destructor
  42 BamMultiReader::~BamMultiReader(void) {
  43     Close();
  44 }
  45
  46 // close the BAM files
  47 void BamMultiReader::Close(void) {
  48
  49     // close all BAM readers and clean up pointers
  50     vector<pair<BamReader*, BamAlignment*> >::iterator readerIter = readers.begin();
  51     vector<pair<BamReader*, BamAlignment*> >::iterator readerEnd  = readers.end();
  52     for ( ; readerIter != readerEnd; ++readerIter) {
  53
  54         BamReader* reader = (*readerIter).first;
  55         BamAlignment* alignment = (*readerIter).second;
  56
  57         // close the reader
  58         if ( reader) reader->Close();
  59
  60         // delete reader pointer
  61         delete reader;
  62         reader = 0;
  63
  64         // delete alignment pointer
  65         delete alignment;
  66         alignment = 0;
  67     }
  68
  69     // clear out the container
  70     readers.clear();
  71 }
  72
  73 // saves index data to BAM index files (".bai"/".bti") where necessary, returns success/fail
  74 bool BamMultiReader::CreateIndexes(bool useStandardIndex) {
  75     bool result = true;
  76     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
  77         BamReader* reader = it->first;
  78         result &= reader->CreateIndex(useStandardIndex);
  79     }
  80     return result;
  81 }
  82
  83 // sets the index caching mode on the readers
  84 void BamMultiReader::SetIndexCacheMode(const BamIndex::BamIndexCacheMode mode) {
  85     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
  86         BamReader* reader = it->first;
  87         reader->SetIndexCacheMode(mode);
  88     }
  89 }
  90
  91 // for debugging
  92 void BamMultiReader::DumpAlignmentIndex(void) {
  93     for (AlignmentIndex::const_iterator it = alignments.begin(); it != alignments.end(); ++it) {
  94         cerr << it->first.first << ":" << it->first.second << " " << it->second.first->GetFilename() << endl;
  95     }
  96 }
  97
  98 // makes a virtual, unified header for all the bam files in the multireader
  99 const string BamMultiReader::GetHeaderText(void) const {
 100
 101     string mergedHeader = "";
 102     map<string, bool> readGroups;
 103
 104     // foreach extraction entry (each BAM file)
 105     for (vector<pair<BamReader*, BamAlignment*> >::const_iterator rs = readers.begin(); rs != readers.end(); ++rs) {
 106
 107         BamReader* reader = rs->first;
 108         string headerText = reader->GetHeaderText();
 109         if ( headerText.empty() ) continue;
 110
 111         map<string, bool> currentFileReadGroups;
 112         stringstream header(headerText);
 113         vector<string> lines;
 114         string item;
 115         while (getline(header, item))
 116             lines.push_back(item);
 117
 118         for (vector<string>::const_iterator it = lines.begin(); it != lines.end(); ++it) {
 119
 120             // get next line from header, skip if empty
 121             string headerLine = *it;
 122             if ( headerLine.empty() ) { continue; }
 123
 124             // if first file, save HD & SQ entries
 125             if ( rs == readers.begin() ) {
 126                 if ( headerLine.find("@HD") == 0 || headerLine.find("@SQ") == 0) {
 127                     mergedHeader.append(headerLine.c_str());
 128                     mergedHeader.append(1, '\n');
 129                 }
 130             }
 131
 132             // (for all files) append RG entries if they are unique
 133             if ( headerLine.find("@RG") == 0 ) {
 134                 stringstream headerLineSs(headerLine);
 135                 string part, readGroupPart, readGroup;
 136                 while(std::getline(headerLineSs, part, '\t')) {
 137                     stringstream partSs(part);
 138                     string subtag;
 139                     std::getline(partSs, subtag, ':');
 140                     if (subtag == "ID") {
 141                         std::getline(partSs, readGroup, ':');
 142                         break;
 143                     }
 144                 }
 145                 if (readGroups.find(readGroup) == readGroups.end()) { // prevents duplicate @RG entries
 146                     mergedHeader.append(headerLine.c_str() );
 147                     mergedHeader.append(1, '\n');
 148                     readGroups[readGroup] = true;
 149                     currentFileReadGroups[readGroup] = true;
 150                 } else {
 151                     // warn iff we are reading one file and discover duplicated @RG tags in the header
 152                     // otherwise, we emit no warning, as we might be merging multiple BAM files with identical @RG tags
 153                     if (currentFileReadGroups.find(readGroup) != currentFileReadGroups.end()) {
 154                         cerr << "WARNING: duplicate @RG tag " << readGroup
 155                             << " entry in header of " << reader->GetFilename() << endl;
 156                     }
 157                 }
 158             }
 159         }
 160     }
 161
 162     // return merged header text
 163     return mergedHeader;
 164 }
 165
 166 // get next alignment among all files
 167 bool BamMultiReader::GetNextAlignment(BamAlignment& nextAlignment) {
 168
 169     // bail out if we are at EOF in all files, means no more alignments to process
 170     if (!HasOpenReaders())
 171         return false;
 172
 173     // when all alignments have stepped into a new target sequence, update our
 174     // current reference sequence id
 175     UpdateReferenceID();
 176
 177     // our lowest alignment and reader will be at the front of our alignment index
 178     BamAlignment* alignment = alignments.begin()->second.second;
 179     BamReader* reader = alignments.begin()->second.first;
 180
 181     // now that we have the lowest alignment in the set, save it by copy to our argument
 182     nextAlignment = BamAlignment(*alignment);
 183
 184     // remove this alignment index entry from our alignment index
 185     alignments.erase(alignments.begin());
 186
 187     // and add another entry if we can get another alignment from the reader
 188     if (reader->GetNextAlignment(*alignment)) {
 189         alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position),
 190                                     make_pair(reader, alignment)));
 191     } else { // do nothing
 192         //cerr << "reached end of file " << lowestReader->GetFilename() << endl;
 193     }
 194
 195     return true;
 196
 197 }
 198
 199 // get next alignment among all files without parsing character data from alignments
 200 bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) {
 201
 202     // bail out if we are at EOF in all files, means no more alignments to process
 203     if (!HasOpenReaders())
 204         return false;
 205
 206     // when all alignments have stepped into a new target sequence, update our
 207     // current reference sequence id
 208     UpdateReferenceID();
 209
 210     // our lowest alignment and reader will be at the front of our alignment index
 211     BamAlignment* alignment = alignments.begin()->second.second;
 212     BamReader* reader = alignments.begin()->second.first;
 213
 214     // now that we have the lowest alignment in the set, save it by copy to our argument
 215     nextAlignment = BamAlignment(*alignment);
 216     //memcpy(&nextAlignment, alignment, sizeof(BamAlignment));
 217
 218     // remove this alignment index entry from our alignment index
 219     alignments.erase(alignments.begin());
 220
 221     // and add another entry if we can get another alignment from the reader
 222     if (reader->GetNextAlignmentCore(*alignment)) {
 223         alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position),
 224                                     make_pair(reader, alignment)));
 225     } else { // do nothing
 226         //cerr << "reached end of file " << lowestReader->GetFilename() << endl;
 227     }
 228
 229     return true;
 230
 231 }
 232
 233 // ---------------------------------------------------------------------------------------
 234 //
 235 // NB: The following GetReferenceX() functions assume that we have identical
 236 // references for all BAM files.  We enforce this by invoking the above
 237 // validation function (ValidateReaders) to verify that our reference data
 238 // is the same across all files on Open, so we will not encounter a situation
 239 // in which there is a mismatch and we are still live.
 240 //
 241 // ---------------------------------------------------------------------------------------
 242
 243 // returns the number of reference sequences
 244 const int BamMultiReader::GetReferenceCount(void) const {
 245     return readers.front().first->GetReferenceCount();
 246 }
 247
 248 // returns vector of reference objects
 249 const BamTools::RefVector BamMultiReader::GetReferenceData(void) const {
 250     return readers.front().first->GetReferenceData();
 251 }
 252
 253 // returns refID from reference name
 254 const int BamMultiReader::GetReferenceID(const string& refName) const {
 255     return readers.front().first->GetReferenceID(refName);
 256 }
 257
 258 // ---------------------------------------------------------------------------------------
 259
 260 // checks if any readers still have alignments
 261 bool BamMultiReader::HasOpenReaders() {
 262     return alignments.size() > 0;
 263 }
 264
 265 // returns whether underlying BAM readers ALL have an index loaded
 266 // this is useful to indicate whether Jump() or SetRegion() are possible
 267 bool BamMultiReader::IsIndexLoaded(void) const {
 268     bool ok = true;
 269     vector<pair<BamReader*, BamAlignment*> >::const_iterator readerIter = readers.begin();
 270     vector<pair<BamReader*, BamAlignment*> >::const_iterator readerEnd  = readers.end();
 271     for ( ; readerIter != readerEnd; ++readerIter ) {
 272         const BamReader* reader = (*readerIter).first;
 273         if ( reader ) ok &= reader->IsIndexLoaded();
 274     }
 275     return ok;
 276 }
 277
 278 // jumps to specified region(refID, leftBound) in BAM files, returns success/fail
 279 bool BamMultiReader::Jump(int refID, int position) {
 280
 281     //if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {
 282     CurrentRefID = refID;
 283     CurrentLeft  = position;
 284
 285     bool result = true;
 286     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
 287         BamReader* reader = it->first;
 288         result &= reader->Jump(refID, position);
 289         if (!result) {
 290             cerr << "ERROR: could not jump " << reader->GetFilename() << " to " << refID << ":" << position << endl;
 291             exit(1);
 292         }
 293     }
 294     if (result) UpdateAlignments();
 295     return result;
 296 }
 297
 298 // opens BAM files
 299 bool BamMultiReader::Open(const vector<string>& filenames, bool openIndexes, bool coreMode, bool preferStandardIndex) {
 300
 301     // for filename in filenames
 302     fileNames = filenames; // save filenames in our multireader
 303     for (vector<string>::const_iterator it = filenames.begin(); it != filenames.end(); ++it) {
 304
 305         const string filename = *it;
 306         BamReader* reader = new BamReader;
 307
 308         bool openedOK = true;
 309         if (openIndexes) {
 310
 311             // leave index filename empty
 312             // this allows BamReader & BamIndex to search for any available
 313             // useStandardIndex gives hint to prefer BAI over BTI
 314             openedOK = reader->Open(filename, "", true, preferStandardIndex);
 315         }
 316
 317         // ignoring index file(s)
 318         else openedOK = reader->Open(filename);
 319
 320         // if file opened ok, check that it can be read
 321         if ( openedOK ) {
 322
 323             bool fileOK = true;
 324             BamAlignment* alignment = new BamAlignment;
 325             fileOK &= ( coreMode ? reader->GetNextAlignmentCore(*alignment) : reader->GetNextAlignment(*alignment) );
 326
 327             if (fileOK) {
 328                 readers.push_back(make_pair(reader, alignment)); // store pointers to our readers for cleanup
 329                 alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position),
 330                                             make_pair(reader, alignment)));
 331             } else {
 332                 cerr << "WARNING: could not read first alignment in " << filename << ", ignoring file" << endl;
 333                 // if only file available & could not be read, return failure
 334                 if ( filenames.size() == 1 ) return false;
 335             }
 336         }
 337
 338         // TODO; any further error handling when openedOK is false ??
 339         else
 340             return false;
 341     }
 342
 343     // files opened ok, at least one alignment could be read,
 344     // now need to check that all files use same reference data
 345     ValidateReaders();
 346     return true;
 347 }
 348
 349 void BamMultiReader::PrintFilenames(void) {
 350     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
 351         BamReader* reader = it->first;
 352         cout << reader->GetFilename() << endl;
 353     }
 354 }
 355
 356 // returns BAM file pointers to beginning of alignment data
 357 bool BamMultiReader::Rewind(void) {
 358     bool result = true;
 359     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
 360         BamReader* reader = it->first;
 361         result &= reader->Rewind();
 362     }
 363     return result;
 364 }
 365
 366 bool BamMultiReader::SetRegion(const int& leftRefID, const int& leftPosition, const int& rightRefID, const int& rightPosition) {
 367     BamRegion region(leftRefID, leftPosition, rightRefID, rightPosition);
 368     return SetRegion(region);
 369 }
 370
 371 bool BamMultiReader::SetRegion(const BamRegion& region) {
 372
 373     Region = region;
 374
 375     // NB: While it may make sense to track readers in which we can
 376     // successfully SetRegion, In practice a failure of SetRegion means "no
 377     // alignments here."  It makes sense to simply accept the failure,
 378     // UpdateAlignments(), and continue.
 379
 380     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
 381         if (!it->first->SetRegion(region)) {
 382             cerr << "ERROR: could not jump " << it->first->GetFilename() << " to "
 383                 << region.LeftRefID << ":" << region.LeftPosition
 384                 << ".." << region.RightRefID << ":" << region.RightPosition << endl;
 385         }
 386     }
 387
 388     UpdateAlignments();
 389     return true;
 390 }
 391
 392 void BamMultiReader::UpdateAlignments(void) {
 393     // Update Alignments
 394     alignments.clear();
 395     for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
 396         BamReader* br = it->first;
 397         BamAlignment* ba = it->second;
 398         if (br->GetNextAlignment(*ba)) {
 399             alignments.insert(make_pair(make_pair(ba->RefID, ba->Position),
 400                                         make_pair(br, ba)));
 401         } else {
 402             // assume BamReader end of region / EOF
 403         }
 404     }
 405 }
 406
 407 // updates the reference id stored in the BamMultiReader
 408 // to reflect the current state of the readers
 409 void BamMultiReader::UpdateReferenceID(void) {
 410     // the alignments are sorted by position, so the first alignment will always have the lowest reference ID
 411     if (alignments.begin()->second.second->RefID != CurrentRefID) {
 412         // get the next reference id
 413         // while there aren't any readers at the next ref id
 414         // increment the ref id
 415         int nextRefID = CurrentRefID;
 416         while (alignments.begin()->second.second->RefID != nextRefID) {
 417             ++nextRefID;
 418         }
 419         //cerr << "updating reference id from " << CurrentRefID << " to " << nextRefID << endl;
 420         CurrentRefID = nextRefID;
 421     }
 422 }
 423
 424 // ValidateReaders checks that all the readers point to BAM files representing
 425 // alignments against the same set of reference sequences, and that the
 426 // sequences are identically ordered.  If these checks fail the operation of
 427 // the multireader is undefined, so we force program exit.
 428 void BamMultiReader::ValidateReaders(void) const {
 429     int firstRefCount = readers.front().first->GetReferenceCount();
 430     BamTools::RefVector firstRefData = readers.front().first->GetReferenceData();
 431     for (vector<pair<BamReader*, BamAlignment*> >::const_iterator it = readers.begin(); it != readers.end(); ++it) {
 432         BamReader* reader = it->first;
 433         BamTools::RefVector currentRefData = reader->GetReferenceData();
 434         BamTools::RefVector::const_iterator f = firstRefData.begin();
 435         BamTools::RefVector::const_iterator c = currentRefData.begin();
 436         if (reader->GetReferenceCount() != firstRefCount || firstRefData.size() != currentRefData.size()) {
 437             cerr << "ERROR: mismatched number of references in " << reader->GetFilename()
 438                       << " expected " << firstRefCount
 439                       << " reference sequences but only found " << reader->GetReferenceCount() << endl;
 440             exit(1);
 441         }
 442         // this will be ok; we just checked above that we have identically-sized sets of references
 443         // here we simply check if they are all, in fact, equal in content
 444         while (f != firstRefData.end()) {
 445             if (f->RefName != c->RefName || f->RefLength != c->RefLength) {
 446                 cerr << "ERROR: mismatched references found in " << reader->GetFilename()
 447                           << " expected: " << endl;
 448                 for (BamTools::RefVector::const_iterator a = firstRefData.begin(); a != firstRefData.end(); ++a)
 449                     cerr << a->RefName << " " << a->RefLength << endl;
 450                 cerr << "but found: " << endl;
 451                 for (BamTools::RefVector::const_iterator a = currentRefData.begin(); a != currentRefData.end(); ++a)
 452                     cerr << a->RefName << " " << a->RefLength << endl;
 453                 exit(1);
 454             }
 455             ++f; ++c;
 456         }
 457     }
 458 }