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