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
10 // ---------------------------------------------------------------------------
11 // Functionality for simultaneously reading multiple BAM files.
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 // ***************************************************************************
29 #include "BamMultiReader.h"
30 using namespace BamTools;
33 // -----------------------------------------------------
34 // BamMultiReader implementation
35 // -----------------------------------------------------
38 BamMultiReader::BamMultiReader(void)
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) {
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
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) {
73 //cerr << "updating reference id from " << CurrentRefID << " to " << nextRefID << endl;
74 CurrentRefID = nextRefID;
78 // checks if any readers still have alignments
79 bool BamMultiReader::HasOpenReaders() {
80 return alignments.size() > 0;
83 // get next alignment among all files
84 bool BamMultiReader::GetNextAlignment(BamAlignment& nextAlignment) {
86 // bail out if we are at EOF in all files, means no more alignments to process
87 if (!HasOpenReaders())
90 // when all alignments have stepped into a new target sequence, update our
91 // current reference sequence id
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;
98 // now that we have the lowest alignment in the set, save it by copy to our argument
99 nextAlignment = BamAlignment(*alignment);
101 // remove this alignment index entry from our alignment index
102 alignments.erase(alignments.begin());
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;
116 // get next alignment among all files without parsing character data from alignments
117 bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) {
119 // bail out if we are at EOF in all files, means no more alignments to process
120 if (!HasOpenReaders())
123 // when all alignments have stepped into a new target sequence, update our
124 // current reference sequence id
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;
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));
135 // remove this alignment index entry from our alignment index
136 alignments.erase(alignments.begin());
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;
150 // jumps to specified region(refID, leftBound) in BAM files, returns success/fail
151 bool BamMultiReader::Jump(int refID, int position) {
153 //if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {
154 CurrentRefID = refID;
155 CurrentLeft = position;
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);
162 cerr << "ERROR: could not jump " << reader->GetFilename() << " to " << refID << ":" << position << endl;
166 if (result) UpdateAlignments();
170 bool BamMultiReader::SetRegion(const int& leftRefID, const int& leftPosition, const int& rightRefID, const int& rightPosition) {
172 BamRegion region(leftRefID, leftPosition, rightRefID, rightPosition);
174 return SetRegion(region);
178 bool BamMultiReader::SetRegion(const BamRegion& region) {
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.
187 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
188 it->first->SetRegion(region);
197 void BamMultiReader::UpdateAlignments(void) {
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),
207 // assume BamReader end of region / EOF
213 bool BamMultiReader::Open(const vector<string> filenames, bool openIndexes, bool coreMode, bool useDefaultIndex) {
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;
221 bool openedOK = true;
224 openedOK = reader->Open(filename, filename + ".bai");
226 openedOK = reader->Open(filename, filename + ".bti");
228 openedOK = reader->Open(filename); // for merging, jumping is disallowed
231 // if file opened ok, check that it can be read
235 BamAlignment* alignment = new BamAlignment;
237 fileOK &= reader->GetNextAlignmentCore(*alignment);
239 fileOK &= reader->GetNextAlignment(*alignment);
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)));
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;
254 // TODO; any more error handling on openedOK ??
259 // files opened ok, at least one alignment could be read,
260 // now need to check that all files use same reference data
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;
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;
279 // returns BAM file pointers to beginning of alignment data
280 bool BamMultiReader::Rewind(void) {
282 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
283 BamReader* reader = it->first;
284 result &= reader->Rewind();
289 // saves index data to BAM index files (".bai") where necessary, returns success/fail
290 bool BamMultiReader::CreateIndexes(bool useDefaultIndex) {
292 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
293 BamReader* reader = it->first;
294 result &= reader->CreateIndex(useDefaultIndex);
299 // makes a virtual, unified header for all the bam files in the multireader
300 const string BamMultiReader::GetHeaderText(void) const {
302 string mergedHeader = "";
303 map<string, bool> readGroups;
305 // foreach extraction entry (each BAM file)
306 for (vector<pair<BamReader*, BamAlignment*> >::const_iterator rs = readers.begin(); rs != readers.end(); ++rs) {
308 map<string, bool> currentFileReadGroups;
310 BamReader* reader = rs->first;
312 stringstream header(reader->GetHeaderText());
313 vector<string> lines;
315 while (getline(header, item))
316 lines.push_back(item);
318 for (vector<string>::const_iterator it = lines.begin(); it != lines.end(); ++it) {
320 // get next line from header, skip if empty
321 string headerLine = *it;
322 if ( headerLine.empty() ) { continue; }
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');
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);
339 std::getline(partSs, subtag, ':');
340 if (subtag == "ID") {
341 std::getline(partSs, readGroup, ':');
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;
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;
362 // return merged header text
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;
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;
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.
408 // returns the number of reference sequences
409 const int BamMultiReader::GetReferenceCount(void) const {
410 return readers.front().first->GetReferenceCount();
413 // returns vector of reference objects
414 const BamTools::RefVector BamMultiReader::GetReferenceData(void) const {
415 return readers.front().first->GetReferenceData();
418 const int BamMultiReader::GetReferenceID(const string& refName) const {
419 return readers.front().first->GetReferenceID(refName);