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