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