1 // ***************************************************************************
2 // bamtools_sort.cpp (c) 2010 Derek Barnett, Erik Garrison
3 // Marth Lab, Department of Biology, Boston College
4 // All rights reserved.
5 // ---------------------------------------------------------------------------
6 // Last modified: 14 June 2011 (DB)
7 // ---------------------------------------------------------------------------
8 // Sorts an input BAM file
9 // ***************************************************************************
11 #include "bamtools_sort.h"
13 #include <api/SamConstants.h>
14 #include <api/BamMultiReader.h>
15 #include <api/BamWriter.h>
16 #include <utils/bamtools_options.h>
17 using namespace BamTools;
31 // ** These defaults should be tweaked & 'optimized' per testing ** //
33 // I say 'optimized' because each system will naturally perform
34 // differently. We will attempt to determine a sensible
35 // compromise that should perform well on average.
36 const unsigned int SORT_DEFAULT_MAX_BUFFER_COUNT = 500000; // max numberOfAlignments for buffer
37 const unsigned int SORT_DEFAULT_MAX_BUFFER_MEMORY = 1024; // Mb
39 // -----------------------------------
40 // comparison objects (for sorting)
42 struct SortLessThanPosition {
43 bool operator() (const BamAlignment& lhs, const BamAlignment& rhs) {
44 if ( lhs.RefID != rhs.RefID )
45 return lhs.RefID < rhs.RefID;
47 return lhs.Position < rhs.Position;
51 struct SortLessThanName {
52 bool operator() (const BamAlignment& lhs, const BamAlignment& rhs) {
53 return lhs.Name < rhs.Name;
57 } // namespace BamTools
59 // ---------------------------------------------
60 // SortSettings implementation
62 struct SortTool::SortSettings {
65 bool HasInputBamFilename;
66 bool HasMaxBufferCount;
67 bool HasMaxBufferMemory;
68 bool HasOutputBamFilename;
72 string InputBamFilename;
73 string OutputBamFilename;
76 unsigned int MaxBufferCount;
77 unsigned int MaxBufferMemory;
81 : HasInputBamFilename(false)
82 , HasMaxBufferCount(false)
83 , HasMaxBufferMemory(false)
84 , HasOutputBamFilename(false)
85 , IsSortingByName(false)
86 , InputBamFilename(Options::StandardIn())
87 , OutputBamFilename(Options::StandardOut())
88 , MaxBufferCount(SORT_DEFAULT_MAX_BUFFER_COUNT)
89 , MaxBufferMemory(SORT_DEFAULT_MAX_BUFFER_MEMORY)
93 // ---------------------------------------------
94 // SortToolPrivate implementation
96 class SortTool::SortToolPrivate {
100 SortToolPrivate(SortTool::SortSettings* settings);
101 ~SortToolPrivate(void) { }
103 // 'public' interface
109 bool CreateSortedTempFile(vector<BamAlignment>& buffer);
110 bool GenerateSortedRuns(void);
111 bool MergeSortedRuns(void);
112 bool WriteTempFile(const vector<BamAlignment>& buffer, const string& tempFilename);
113 void SortBuffer(vector<BamAlignment>& buffer);
117 SortTool::SortSettings* m_settings;
118 string m_tempFilenameStub;
121 RefVector m_references;
122 vector<string> m_tempFilenames;
126 SortTool::SortToolPrivate::SortToolPrivate(SortTool::SortSettings* settings)
127 : m_settings(settings)
130 // set filename stub depending on inputfile path
131 // that way multiple sort runs don't trip on each other's temp files
133 size_t extensionFound = m_settings->InputBamFilename.find(".bam");
134 if (extensionFound != string::npos )
135 m_tempFilenameStub = m_settings->InputBamFilename.substr(0,extensionFound);
136 m_tempFilenameStub.append(".sort.temp.");
140 // generates mutiple sorted temp BAM files from single unsorted BAM file
141 bool SortTool::SortToolPrivate::GenerateSortedRuns(void) {
143 // open input BAM file
145 if ( !reader.Open(m_settings->InputBamFilename) ) {
146 cerr << "bamtools sort ERROR: could not open " << m_settings->InputBamFilename
147 << " for reading... Aborting." << endl;
151 // get basic data that will be shared by all temp/output files
152 SamHeader header = reader.GetHeader();
153 header.SortOrder = ( m_settings->IsSortingByName
154 ? Constants::SAM_HD_SORTORDER_QUERYNAME
155 : Constants::SAM_HD_SORTORDER_COORDINATE );
156 m_headerText = header.ToString();
157 m_references = reader.GetReferenceData();
159 // set up alignments buffer
161 vector<BamAlignment> buffer;
162 buffer.reserve( (size_t)(m_settings->MaxBufferCount*1.1) );
163 bool bufferFull = false;
165 // if sorting by name, we need to generate full char data
166 // so can't use GetNextAlignmentCore()
167 if ( m_settings->IsSortingByName ) {
169 // iterate through file
170 while ( reader.GetNextAlignment(al)) {
172 // check buffer's usage
173 bufferFull = ( buffer.size() >= m_settings->MaxBufferCount );
175 // store alignments until buffer is "full"
177 buffer.push_back(al);
179 // if buffer is "full"
182 // push any unmapped reads into buffer,
183 // don't want to split these into a separate temp file
184 if ( !al.IsMapped() )
185 buffer.push_back(al);
187 // "al" is mapped, so create a sorted temp file with current buffer contents
188 // then push "al" into fresh buffer
190 CreateSortedTempFile(buffer);
191 buffer.push_back(al);
197 // sorting by position, can take advantage of GNACore() speedup
200 // iterate through file
201 while ( reader.GetNextAlignmentCore(al) ) {
203 // check buffer's usage
204 bufferFull = ( buffer.size() >= m_settings->MaxBufferCount );
206 // store alignments until buffer is "full"
208 buffer.push_back(al);
210 // if buffer is "full"
213 // push any unmapped reads into buffer,
214 // don't want to split these into a separate temp file
215 if ( !al.IsMapped() )
216 buffer.push_back(al);
218 // "al" is mapped, so create a sorted temp file with current buffer contents
219 // then push "al" into fresh buffer
221 CreateSortedTempFile(buffer);
222 buffer.push_back(al);
228 // handle any leftover buffer contents
229 if ( !buffer.empty() )
230 CreateSortedTempFile(buffer);
232 // close reader & return success
237 bool SortTool::SortToolPrivate::CreateSortedTempFile(vector<BamAlignment>& buffer) {
242 // write sorted contents to temp file, store success/fail
243 stringstream tempStr;
244 tempStr << m_tempFilenameStub << m_numberOfRuns;
245 bool success = WriteTempFile( buffer, tempStr.str() );
247 // save temp filename for merging later
248 m_tempFilenames.push_back(tempStr.str());
250 // clear buffer contents & update run counter
254 // return success/fail of writing to temp file
255 // TODO: a failure returned here is not actually caught and handled anywhere
259 // merges sorted temp BAM files into single sorted output BAM file
260 bool SortTool::SortToolPrivate::MergeSortedRuns(void) {
262 // open up multi reader for all of our temp files
263 // this might get broken up if we do a multi-pass system later ??
264 BamMultiReader multiReader;
265 if ( !multiReader.Open(m_tempFilenames) ) {
266 cerr << "bamtools sort ERROR: could not open BamMultiReader for merging temp files... Aborting."
271 // set sort order for merge
272 if ( m_settings->IsSortingByName )
273 multiReader.SetSortOrder(BamMultiReader::SortedByReadName);
275 multiReader.SetSortOrder(BamMultiReader::SortedByPosition);
277 // open writer for our completely sorted output BAM file
278 BamWriter mergedWriter;
279 if ( !mergedWriter.Open(m_settings->OutputBamFilename, m_headerText, m_references) ) {
280 cerr << "bamtools sort ERROR: could not open " << m_settings->OutputBamFilename
281 << " for writing... Aborting." << endl;
286 // while data available in temp files
288 while ( multiReader.GetNextAlignmentCore(al) )
289 mergedWriter.SaveAlignment(al);
293 mergedWriter.Close();
295 // delete all temp files
296 vector<string>::const_iterator tempIter = m_tempFilenames.begin();
297 vector<string>::const_iterator tempEnd = m_tempFilenames.end();
298 for ( ; tempIter != tempEnd; ++tempIter ) {
299 const string& tempFilename = (*tempIter);
300 remove(tempFilename.c_str());
306 bool SortTool::SortToolPrivate::Run(void) {
308 // this does a single pass, chunking up the input file into smaller sorted temp files,
309 // then write out using BamMultiReader to handle merging
311 if ( GenerateSortedRuns() )
312 return MergeSortedRuns();
317 void SortTool::SortToolPrivate::SortBuffer(vector<BamAlignment>& buffer) {
319 // ** add further custom sort options later ?? **
321 // sort buffer by desired method
322 if ( m_settings->IsSortingByName )
323 std::stable_sort ( buffer.begin(), buffer.end(), SortLessThanName() );
325 std::stable_sort ( buffer.begin(), buffer.end(), SortLessThanPosition() );
328 bool SortTool::SortToolPrivate::WriteTempFile(const vector<BamAlignment>& buffer,
329 const string& tempFilename)
331 // open temp file for writing
332 BamWriter tempWriter;
333 if ( !tempWriter.Open(tempFilename, m_headerText, m_references) ) {
334 cerr << "bamtools sort ERROR: could not open " << tempFilename
335 << " for writing." << endl;
340 vector<BamAlignment>::const_iterator buffIter = buffer.begin();
341 vector<BamAlignment>::const_iterator buffEnd = buffer.end();
342 for ( ; buffIter != buffEnd; ++buffIter ) {
343 const BamAlignment& al = (*buffIter);
344 tempWriter.SaveAlignment(al);
347 // close temp file & return success
352 // ---------------------------------------------
353 // SortTool implementation
355 SortTool::SortTool(void)
357 , m_settings(new SortSettings)
360 // set program details
361 Options::SetProgramInfo("bamtools sort", "sorts a BAM file", "[-in <filename>] [-out <filename>] [sortOptions]");
364 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
365 Options::AddValueOption("-in", "BAM filename", "the input BAM file", "",
366 m_settings->HasInputBamFilename, m_settings->InputBamFilename,
367 IO_Opts, Options::StandardIn());
368 Options::AddValueOption("-out", "BAM filename", "the output BAM file", "",
369 m_settings->HasOutputBamFilename, m_settings->OutputBamFilename,
370 IO_Opts, Options::StandardOut());
372 OptionGroup* SortOpts = Options::CreateOptionGroup("Sorting Methods");
373 Options::AddOption("-byname", "sort by alignment name", m_settings->IsSortingByName, SortOpts);
375 OptionGroup* MemOpts = Options::CreateOptionGroup("Memory Settings");
376 Options::AddValueOption("-n", "count", "max number of alignments per tempfile", "",
377 m_settings->HasMaxBufferCount, m_settings->MaxBufferCount,
378 MemOpts, SORT_DEFAULT_MAX_BUFFER_COUNT);
379 Options::AddValueOption("-mem", "Mb", "max memory to use", "",
380 m_settings->HasMaxBufferMemory, m_settings->MaxBufferMemory,
381 MemOpts, SORT_DEFAULT_MAX_BUFFER_MEMORY);
384 SortTool::~SortTool(void) {
393 int SortTool::Help(void) {
394 Options::DisplayHelp();
398 int SortTool::Run(int argc, char* argv[]) {
400 // parse command line arguments
401 Options::Parse(argc, argv, 1);
403 // initialize SortTool with settings
404 m_impl = new SortToolPrivate(m_settings);
406 // run SortTool, return success/fail