src/toolkit/bamtools_split.cpp

   1 // ***************************************************************************
   2 // bamtools_split.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 7 April 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Splits a BAM file on user-specified property, creating a new BAM output
   8 // file for each value found
   9 // ***************************************************************************
  10
  11 #include "bamtools_split.h"
  12
  13 #include <api/BamConstants.h>
  14 #include <api/BamReader.h>
  15 #include <api/BamWriter.h>
  16 #include <utils/bamtools_options.h>
  17 #include <utils/bamtools_variant.h>
  18 using namespace BamTools;
  19
  20 #include <ctime>
  21 #include <iostream>
  22 #include <map>
  23 #include <sstream>
  24 #include <string>
  25 #include <vector>
  26 using namespace std;
  27
  28 namespace BamTools {
  29
  30 // string constants
  31 static const string SPLIT_MAPPED_TOKEN    = ".MAPPED";
  32 static const string SPLIT_UNMAPPED_TOKEN  = ".UNMAPPED";
  33 static const string SPLIT_PAIRED_TOKEN    = ".PAIRED_END";
  34 static const string SPLIT_SINGLE_TOKEN    = ".SINGLE_END";
  35 static const string SPLIT_REFERENCE_TOKEN = ".REF_";
  36
  37 string GetTimestampString(void) {
  38
  39     // get human readable timestamp
  40     time_t currentTime;
  41     time(&currentTime);
  42     stringstream timeStream("");
  43     timeStream << ctime(&currentTime);
  44
  45     // convert whitespace to '_'
  46     string timeString = timeStream.str();
  47     size_t found = timeString.find(" ");
  48     while (found != string::npos) {
  49         timeString.replace(found, 1, "_");
  50         found = timeString.find(" ", found+1);
  51     }
  52     return timeString;
  53 }
  54
  55 // remove copy of filename without extension
  56 // (so /path/to/file.txt becomes /path/to/file )
  57 string RemoveFilenameExtension(const string& filename) {
  58     size_t found = filename.rfind(".");
  59     return filename.substr(0, found);
  60 }
  61
  62 } // namespace BamTools
  63
  64 // ---------------------------------------------
  65 // SplitSettings implementation
  66
  67 struct SplitTool::SplitSettings {
  68
  69     // flags
  70     bool HasInputFilename;
  71     bool HasCustomOutputStub;
  72     bool IsSplittingMapped;
  73     bool IsSplittingPaired;
  74     bool IsSplittingReference;
  75     bool IsSplittingTag;
  76
  77     // string args
  78     string CustomOutputStub;
  79     string InputFilename;
  80     string TagToSplit;
  81
  82     // constructor
  83     SplitSettings(void)
  84         : HasInputFilename(false)
  85         , HasCustomOutputStub(false)
  86         , IsSplittingMapped(false)
  87         , IsSplittingPaired(false)
  88         , IsSplittingReference(false)
  89         , IsSplittingTag(false)
  90         , CustomOutputStub("")
  91         , InputFilename(Options::StandardIn())
  92         , TagToSplit("")
  93     { }
  94 };
  95
  96 // ---------------------------------------------
  97 // SplitToolPrivate declaration
  98
  99 class SplitTool::SplitToolPrivate {
 100
 101     // ctor & dtor
 102     public:
 103         SplitToolPrivate(SplitTool::SplitSettings* settings)
 104             : m_settings(settings)
 105         { }
 106
 107         ~SplitToolPrivate(void) {
 108             m_reader.Close();
 109         }
 110
 111     // 'public' interface
 112     public:
 113         bool Run(void);
 114
 115     // internal methods
 116     private:
 117         // close & delete BamWriters in map
 118         template<typename T>
 119         void CloseWriters(map<T, BamWriter*>& writers);
 120         // calculate output stub based on IO args given
 121         void DetermineOutputFilenameStub(void);
 122         // open our BamReader
 123         bool OpenReader(void);
 124         // split alignments in BAM file based on isMapped property
 125         bool SplitMapped(void);
 126         // split alignments in BAM file based on isPaired property
 127         bool SplitPaired(void);
 128         // split alignments in BAM file based on refID property
 129         bool SplitReference(void);
 130         // finds first alignment and calls corresponding SplitTagImpl<>
 131         // depending on tag type
 132         bool SplitTag(void);
 133         // templated split tag implementation
 134         // handle the various types that are possible for tags
 135         template<typename T>
 136         bool SplitTagImpl(BamAlignment& al);
 137
 138     // data members
 139     private:
 140         SplitTool::SplitSettings* m_settings;
 141         string m_outputFilenameStub;
 142         BamReader m_reader;
 143         string m_header;
 144         RefVector m_references;
 145 };
 146
 147 void SplitTool::SplitToolPrivate::DetermineOutputFilenameStub(void) {
 148
 149     // if user supplied output filename stub, use that
 150     if ( m_settings->HasCustomOutputStub )
 151         m_outputFilenameStub = m_settings->CustomOutputStub;
 152
 153     // else if user supplied input BAM filename, use that (minus ".bam" extension) as stub
 154     else if ( m_settings->HasInputFilename )
 155         m_outputFilenameStub = RemoveFilenameExtension(m_settings->InputFilename);
 156
 157     // otherwise, user did not specify -stub, and input is coming from STDIN
 158     // generate stub from timestamp
 159     else m_outputFilenameStub = GetTimestampString();
 160 }
 161
 162 bool SplitTool::SplitToolPrivate::OpenReader(void) {
 163
 164     // attempt to open BAM file
 165     if ( !m_reader.Open(m_settings->InputFilename) ) {
 166         cerr << "bamtools split ERROR: could not open BAM file: " << m_settings->InputFilename << endl;
 167         return false;
 168     }
 169
 170     // save file 'metadata' & return success
 171     m_header     = m_reader.GetHeaderText();
 172     m_references = m_reader.GetReferenceData();
 173     return true;
 174 }
 175
 176 bool SplitTool::SplitToolPrivate::Run(void) {
 177
 178     // determine output stub
 179     DetermineOutputFilenameStub();
 180
 181     // open up BamReader
 182     if ( !OpenReader() )
 183         return false;
 184
 185     // determine split type from settings
 186     if ( m_settings->IsSplittingMapped )    return SplitMapped();
 187     if ( m_settings->IsSplittingPaired )    return SplitPaired();
 188     if ( m_settings->IsSplittingReference ) return SplitReference();
 189     if ( m_settings->IsSplittingTag )       return SplitTag();
 190
 191     // if we get here, no property was specified
 192     cerr << "bamtools split ERROR: no property given to split on... " << endl
 193          << "Please use -mapped, -paired, -reference, or -tag TAG to specifiy desired split behavior." << endl;
 194     return false;
 195 }
 196
 197 bool SplitTool::SplitToolPrivate::SplitMapped(void) {
 198
 199     // set up splitting data structure
 200     map<bool, BamWriter*> outputFiles;
 201     map<bool, BamWriter*>::iterator writerIter;
 202
 203     // iterate through alignments
 204     BamAlignment al;
 205     BamWriter* writer;
 206     bool isCurrentAlignmentMapped;
 207     while ( m_reader.GetNextAlignment(al) ) {
 208
 209         // see if bool value exists
 210         isCurrentAlignmentMapped = al.IsMapped();
 211         writerIter = outputFiles.find(isCurrentAlignmentMapped);
 212
 213         // if no writer associated with this value
 214         if ( writerIter == outputFiles.end() ) {
 215
 216             // open new BamWriter
 217             const string outputFilename = m_outputFilenameStub + ( isCurrentAlignmentMapped
 218                                                                   ? SPLIT_MAPPED_TOKEN
 219                                                                   : SPLIT_UNMAPPED_TOKEN ) + ".bam";
 220             writer = new BamWriter;
 221             if ( !writer->Open(outputFilename, m_header, m_references) ) {
 222                 cerr << "bamtools split ERROR: could not open " << outputFilename
 223                      << " for writing." << endl;
 224                 return false;
 225             }
 226
 227             // store in map
 228             outputFiles.insert( make_pair(isCurrentAlignmentMapped, writer) );
 229         }
 230
 231         // else grab corresponding writer
 232         else writer = (*writerIter).second;
 233
 234         // store alignment in proper BAM output file
 235         if ( writer )
 236             writer->SaveAlignment(al);
 237     }
 238
 239     // clean up BamWriters
 240     CloseWriters(outputFiles);
 241
 242     // return success
 243     return true;
 244 }
 245
 246 bool SplitTool::SplitToolPrivate::SplitPaired(void) {
 247
 248     // set up splitting data structure
 249     map<bool, BamWriter*> outputFiles;
 250     map<bool, BamWriter*>::iterator writerIter;
 251
 252     // iterate through alignments
 253     BamAlignment al;
 254     BamWriter* writer;
 255     bool isCurrentAlignmentPaired;
 256     while ( m_reader.GetNextAlignment(al) ) {
 257
 258         // see if bool value exists
 259         isCurrentAlignmentPaired = al.IsPaired();
 260         writerIter = outputFiles.find(isCurrentAlignmentPaired);
 261
 262         // if no writer associated with this value
 263         if ( writerIter == outputFiles.end() ) {
 264
 265             // open new BamWriter
 266             const string outputFilename = m_outputFilenameStub + ( isCurrentAlignmentPaired
 267                                                                   ? SPLIT_PAIRED_TOKEN
 268                                                                   : SPLIT_SINGLE_TOKEN ) + ".bam";
 269             writer = new BamWriter;
 270             if ( !writer->Open(outputFilename, m_header, m_references) ) {
 271                 cerr << "bamtool split ERROR: could not open " << outputFilename
 272                      << " for writing." << endl;
 273                 return false;
 274             }
 275
 276             // store in map
 277             outputFiles.insert( make_pair(isCurrentAlignmentPaired, writer) );
 278         }
 279
 280         // else grab corresponding writer
 281         else writer = (*writerIter).second;
 282
 283         // store alignment in proper BAM output file
 284         if ( writer )
 285             writer->SaveAlignment(al);
 286     }
 287
 288     // clean up BamWriters
 289     CloseWriters(outputFiles);
 290
 291     // return success
 292     return true;
 293 }
 294
 295 bool SplitTool::SplitToolPrivate::SplitReference(void) {
 296
 297     // set up splitting data structure
 298     map<int32_t, BamWriter*> outputFiles;
 299     map<int32_t, BamWriter*>::iterator writerIter;
 300
 301     // iterate through alignments
 302     BamAlignment al;
 303     BamWriter* writer;
 304     int32_t currentRefId;
 305     while ( m_reader.GetNextAlignment(al) ) {
 306
 307         // see if bool value exists
 308         currentRefId = al.RefID;
 309         writerIter = outputFiles.find(currentRefId);
 310
 311         // if no writer associated with this value
 312         if ( writerIter == outputFiles.end() ) {
 313
 314             // fetch reference name for ID
 315             string refName;
 316             if ( currentRefId == -1 )
 317                 refName = "unmapped";
 318             else
 319                 refName = m_references.at(currentRefId).RefName;
 320
 321             // construct new output filename
 322             const string outputFilename = m_outputFilenameStub + SPLIT_REFERENCE_TOKEN + refName + ".bam";
 323
 324             // open new BamWriter
 325             writer = new BamWriter;
 326             if ( !writer->Open(outputFilename, m_header, m_references) ) {
 327                 cerr << "bamtools split ERROR: could not open " << outputFilename
 328                      << " for writing." << endl;
 329                 return false;
 330             }
 331
 332             // store in map
 333             outputFiles.insert( make_pair(currentRefId, writer) );
 334         }
 335
 336         // else grab corresponding writer
 337         else writer = (*writerIter).second;
 338
 339         // store alignment in proper BAM output file
 340         if ( writer )
 341             writer->SaveAlignment(al);
 342     }
 343
 344     // clean up BamWriters
 345     CloseWriters(outputFiles);
 346
 347     // return success
 348     return true;
 349 }
 350
 351 // finds first alignment and calls corresponding SplitTagImpl<>() depending on tag type
 352 bool SplitTool::SplitToolPrivate::SplitTag(void) {
 353
 354     // iterate through alignments, until we hit TAG
 355     BamAlignment al;
 356     while ( m_reader.GetNextAlignment(al) ) {
 357
 358         // look for tag in this alignment and get tag type
 359         char tagType(0);
 360         if ( !al.GetTagType(m_settings->TagToSplit, tagType) )
 361             continue;
 362
 363         // request split method based on tag type
 364         // pass it the current alignment found
 365         switch ( tagType ) {
 366
 367             case (Constants::BAM_TAG_TYPE_INT8)  :
 368             case (Constants::BAM_TAG_TYPE_INT16) :
 369             case (Constants::BAM_TAG_TYPE_INT32) :
 370                 return SplitTagImpl<int32_t>(al);
 371
 372             case (Constants::BAM_TAG_TYPE_UINT8)  :
 373             case (Constants::BAM_TAG_TYPE_UINT16) :
 374             case (Constants::BAM_TAG_TYPE_UINT32) :
 375                 return SplitTagImpl<uint32_t>(al);
 376
 377             case (Constants::BAM_TAG_TYPE_FLOAT)  :
 378                 return SplitTagImpl<float>(al);
 379
 380             case (Constants::BAM_TAG_TYPE_ASCII)  :
 381             case (Constants::BAM_TAG_TYPE_STRING) :
 382             case (Constants::BAM_TAG_TYPE_HEX)    :
 383                 return SplitTagImpl<string>(al);
 384
 385             case (Constants::BAM_TAG_TYPE_ARRAY) :
 386                 cerr << "bamtools split ERROR: array tag types are not supported" << endl;
 387                 return false;
 388
 389             default:
 390                 cerr << "bamtools split ERROR: unknown tag type encountered: " << tagType << endl;
 391                 return false;
 392         }
 393     }
 394
 395     // tag not found, but that's not an error - return success
 396     return true;
 397 }
 398
 399 // --------------------------------------------------------------------------------
 400 // template method implementation
 401 // *Technical Note* - use of template methods declared & defined in ".cpp" file
 402 //                    goes against normal practices, but works here because these
 403 //                    are purely internal (no one can call from outside this file)
 404
 405 // close BamWriters & delete pointers
 406 template<typename T>
 407 void SplitTool::SplitToolPrivate::CloseWriters(map<T, BamWriter*>& writers) {
 408
 409     typedef map<T, BamWriter*> WriterMap;
 410     typedef typename WriterMap::iterator WriterMapIterator;
 411
 412     // iterate over writers
 413     WriterMapIterator writerIter = writers.begin();
 414     WriterMapIterator writerEnd  = writers.end();
 415     for ( ; writerIter != writerEnd; ++writerIter ) {
 416         BamWriter* writer = (*writerIter).second;
 417         if ( writer == 0 ) continue;
 418
 419         // close BamWriter
 420         writer->Close();
 421
 422         // destroy BamWriter
 423         delete writer;
 424         writer = 0;
 425     }
 426
 427     // clear the container (destroying the items doesn't remove them)
 428     writers.clear();
 429 }
 430
 431 // handle the various types that are possible for tags
 432 template<typename T>
 433 bool SplitTool::SplitToolPrivate::SplitTagImpl(BamAlignment& al) {
 434
 435     typedef T TagValueType;
 436     typedef map<TagValueType, BamWriter*> WriterMap;
 437     typedef typename WriterMap::iterator WriterMapIterator;
 438
 439     // set up splitting data structure
 440     WriterMap outputFiles;
 441     WriterMapIterator writerIter;
 442
 443     // local variables
 444     const string tag = m_settings->TagToSplit;
 445     BamWriter* writer;
 446     stringstream outputFilenameStream("");
 447     TagValueType currentValue;
 448
 449     // retrieve first alignment tag value
 450     if ( al.GetTag(tag, currentValue) ) {
 451
 452         // open new BamWriter, save first alignment
 453         outputFilenameStream << m_outputFilenameStub << ".TAG_" << tag << "_" << currentValue << ".bam";
 454         writer = new BamWriter;
 455         if ( !writer->Open(outputFilenameStream.str(), m_header, m_references) ) {
 456             cerr << "bamtools split ERROR: could not open " << outputFilenameStream.str()
 457                  << " for writing." << endl;
 458             return false;
 459         }
 460         writer->SaveAlignment(al);
 461
 462         // store in map
 463         outputFiles.insert( make_pair(currentValue, writer) );
 464
 465         // reset stream
 466         outputFilenameStream.str("");
 467     }
 468
 469     // iterate through remaining alignments
 470     while ( m_reader.GetNextAlignment(al) ) {
 471
 472         // skip if this alignment doesn't have TAG
 473         if ( !al.GetTag(tag, currentValue) ) continue;
 474
 475         // look up tag value in map
 476         writerIter = outputFiles.find(currentValue);
 477
 478         // if no writer associated with this value
 479         if ( writerIter == outputFiles.end() ) {
 480
 481             // open new BamWriter
 482             outputFilenameStream << m_outputFilenameStub << ".TAG_" << tag << "_" << currentValue << ".bam";
 483             writer = new BamWriter;
 484             if ( !writer->Open(outputFilenameStream.str(), m_header, m_references) ) {
 485                 cerr << "bamtool split ERROR: could not open " << outputFilenameStream.str()
 486                      << " for writing." << endl;
 487                 return false;
 488             }
 489
 490             // store in map
 491             outputFiles.insert( make_pair(currentValue, writer) );
 492
 493             // reset stream
 494             outputFilenameStream.str("");
 495         }
 496
 497         // else grab corresponding writer
 498         else writer = (*writerIter).second;
 499
 500         // store alignment in proper BAM output file
 501         if ( writer )
 502             writer->SaveAlignment(al);
 503     }
 504
 505     // clean up BamWriters
 506     CloseWriters(outputFiles);
 507
 508     // return success
 509     return true;
 510 }
 511
 512 // ---------------------------------------------
 513 // SplitTool implementation
 514
 515 SplitTool::SplitTool(void)
 516     : AbstractTool()
 517     , m_settings(new SplitSettings)
 518     , m_impl(0)
 519 {
 520     // set program details
 521     Options::SetProgramInfo("bamtools split", "splits a BAM file on user-specified property, creating a new BAM output file for each value found", "[-in <filename>] [-stub <filename stub>] < -mapped | -paired | -reference | -tag <TAG> > ");
 522
 523     // set up options
 524     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 525     Options::AddValueOption("-in",   "BAM filename",  "the input BAM file",  "", m_settings->HasInputFilename,  m_settings->InputFilename,  IO_Opts, Options::StandardIn());
 526     Options::AddValueOption("-stub", "filename stub", "prefix stub for output BAM files (default behavior is to use input filename, without .bam extension, as stub). If input is stdin and no stub provided, a timestamp is generated as the stub.", "", m_settings->HasCustomOutputStub, m_settings->CustomOutputStub, IO_Opts);
 527
 528     OptionGroup* SplitOpts = Options::CreateOptionGroup("Split Options");
 529     Options::AddOption("-mapped",    "split mapped/unmapped alignments",       m_settings->IsSplittingMapped,    SplitOpts);
 530     Options::AddOption("-paired",    "split single-end/paired-end alignments", m_settings->IsSplittingPaired,    SplitOpts);
 531     Options::AddOption("-reference", "split alignments by reference",          m_settings->IsSplittingReference, SplitOpts);
 532     Options::AddValueOption("-tag", "tag name", "splits alignments based on all values of TAG encountered (i.e. -tag RG creates a BAM file for each read group in original BAM file)", "",
 533                             m_settings->IsSplittingTag, m_settings->TagToSplit, SplitOpts);
 534 }
 535
 536 SplitTool::~SplitTool(void) {
 537
 538     delete m_settings;
 539     m_settings = 0;
 540
 541     delete m_impl;
 542     m_impl = 0;
 543 }
 544
 545 int SplitTool::Help(void) {
 546     Options::DisplayHelp();
 547     return 0;
 548 }
 549
 550 int SplitTool::Run(int argc, char* argv[]) {
 551
 552     // parse command line arguments
 553     Options::Parse(argc, argv, 1);
 554
 555     // initialize SplitTool with settings
 556     m_impl = new SplitToolPrivate(m_settings);
 557
 558     // run SplitTool, return success/fail
 559     if ( m_impl->Run() )
 560         return 0;
 561     else
 562         return 1;
 563 }