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