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
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 // ***************************************************************************
30 #include "BamMultiReader.h"
31 using namespace BamTools;
36 bool FileExists(const string& filename) {
37 ifstream f(filename.c_str(), ifstream::in);
41 } // namespace BamTools
43 // -----------------------------------------------------
44 // BamMultiReader implementation
45 // -----------------------------------------------------
48 BamMultiReader::BamMultiReader(void)
54 BamMultiReader::~BamMultiReader(void) {
58 // close the BAM files
59 void BamMultiReader::Close(void) {
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) {
66 BamReader* reader = (*readerIter).first;
67 BamAlignment* alignment = (*readerIter).second;
70 if ( reader) reader->Close();
72 // delete reader pointer
76 // delete alignment pointer
81 // clear out the container
85 // saves index data to BAM index files (".bai"/".bti") where necessary, returns success/fail
86 bool BamMultiReader::CreateIndexes(bool useDefaultIndex) {
88 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
89 BamReader* reader = it->first;
90 result &= reader->CreateIndex(useDefaultIndex);
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;
102 // makes a virtual, unified header for all the bam files in the multireader
103 const string BamMultiReader::GetHeaderText(void) const {
105 string mergedHeader = "";
106 map<string, bool> readGroups;
108 // foreach extraction entry (each BAM file)
109 for (vector<pair<BamReader*, BamAlignment*> >::const_iterator rs = readers.begin(); rs != readers.end(); ++rs) {
111 map<string, bool> currentFileReadGroups;
113 BamReader* reader = rs->first;
115 stringstream header(reader->GetHeaderText());
116 vector<string> lines;
118 while (getline(header, item))
119 lines.push_back(item);
121 for (vector<string>::const_iterator it = lines.begin(); it != lines.end(); ++it) {
123 // get next line from header, skip if empty
124 string headerLine = *it;
125 if ( headerLine.empty() ) { continue; }
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');
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);
142 std::getline(partSs, subtag, ':');
143 if (subtag == "ID") {
144 std::getline(partSs, readGroup, ':');
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;
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;
165 // return merged header text
169 // get next alignment among all files
170 bool BamMultiReader::GetNextAlignment(BamAlignment& nextAlignment) {
172 // bail out if we are at EOF in all files, means no more alignments to process
173 if (!HasOpenReaders())
176 // when all alignments have stepped into a new target sequence, update our
177 // current reference sequence id
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;
184 // now that we have the lowest alignment in the set, save it by copy to our argument
185 nextAlignment = BamAlignment(*alignment);
187 // remove this alignment index entry from our alignment index
188 alignments.erase(alignments.begin());
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;
202 // get next alignment among all files without parsing character data from alignments
203 bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) {
205 // bail out if we are at EOF in all files, means no more alignments to process
206 if (!HasOpenReaders())
209 // when all alignments have stepped into a new target sequence, update our
210 // current reference sequence id
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;
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));
221 // remove this alignment index entry from our alignment index
222 alignments.erase(alignments.begin());
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;
236 // ---------------------------------------------------------------------------------------
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.
244 // ---------------------------------------------------------------------------------------
246 // returns the number of reference sequences
247 const int BamMultiReader::GetReferenceCount(void) const {
248 return readers.front().first->GetReferenceCount();
251 // returns vector of reference objects
252 const BamTools::RefVector BamMultiReader::GetReferenceData(void) const {
253 return readers.front().first->GetReferenceData();
256 // returns refID from reference name
257 const int BamMultiReader::GetReferenceID(const string& refName) const {
258 return readers.front().first->GetReferenceID(refName);
261 // ---------------------------------------------------------------------------------------
263 // checks if any readers still have alignments
264 bool BamMultiReader::HasOpenReaders() {
265 return alignments.size() > 0;
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 {
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();
281 // jumps to specified region(refID, leftBound) in BAM files, returns success/fail
282 bool BamMultiReader::Jump(int refID, int position) {
284 //if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {
285 CurrentRefID = refID;
286 CurrentLeft = position;
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);
293 cerr << "ERROR: could not jump " << reader->GetFilename() << " to " << refID << ":" << position << endl;
297 if (result) UpdateAlignments();
302 bool BamMultiReader::Open(const vector<string> filenames, bool openIndexes, bool coreMode, bool useDefaultIndex) {
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) {
308 const string filename = *it;
309 BamReader* reader = new BamReader;
311 bool openedOK = true;
314 const string defaultIndexFilename = filename + ".bai";
315 const string bamToolsIndexFilename = filename + ".bti";
317 // if default BAM index requested and one exists
318 if ( useDefaultIndex && FileExists(defaultIndexFilename) )
319 openedOK = reader->Open(filename, defaultIndexFilename);
321 // else see if BamTools index exists
322 else if ( FileExists(bamToolsIndexFilename) )
323 openedOK = reader->Open(filename, bamToolsIndexFilename);
325 // else none exist... just open without
327 openedOK = reader->Open(filename);
330 // ignoring index file(s)
331 else openedOK = reader->Open(filename);
333 // if file opened ok, check that it can be read
337 BamAlignment* alignment = new BamAlignment;
338 fileOK &= ( coreMode ? reader->GetNextAlignmentCore(*alignment) : reader->GetNextAlignment(*alignment) );
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)));
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;
351 // TODO; any more error handling when openedOKvis false ??
356 // files opened ok, at least one alignment could be read,
357 // now need to check that all files use same reference data
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;
369 // returns BAM file pointers to beginning of alignment data
370 bool BamMultiReader::Rewind(void) {
372 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
373 BamReader* reader = it->first;
374 result &= reader->Rewind();
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);
384 bool BamMultiReader::SetRegion(const BamRegion& region) {
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.
393 for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
394 it->first->SetRegion(region);
401 void BamMultiReader::UpdateAlignments(void) {
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),
411 // assume BamReader end of region / EOF
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) {
428 //cerr << "updating reference id from " << CurrentRefID << " to " << nextRefID << endl;
429 CurrentRefID = nextRefID;
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;
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;