src/api/BamMultiReader.cpp

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