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             // open new BamWriter
 315             const string refName = m_references.at(currentRefId).RefName;
 316             const string outputFilename = m_outputFilenameStub + SPLIT_REFERENCE_TOKEN + refName + ".bam";
 317             writer = new BamWriter;
 318             if ( !writer->Open(outputFilename, m_header, m_references) ) {
 319                 cerr << "bamtools split ERROR: could not open " << outputFilename
 320                      << " for writing." << endl;
 321                 return false;
 322             }
 323
 324             // store in map
 325             outputFiles.insert( make_pair(currentRefId, writer) );
 326         }
 327
 328         // else grab corresponding writer
 329         else writer = (*writerIter).second;
 330
 331         // store alignment in proper BAM output file
 332         if ( writer )
 333             writer->SaveAlignment(al);
 334     }
 335
 336     // clean up BamWriters
 337     CloseWriters(outputFiles);
 338
 339     // return success
 340     return true;
 341 }
 342
 343 // finds first alignment and calls corresponding SplitTagImpl<>() depending on tag type
 344 bool SplitTool::SplitToolPrivate::SplitTag(void) {
 345
 346     // iterate through alignments, until we hit TAG
 347     BamAlignment al;
 348     while ( m_reader.GetNextAlignment(al) ) {
 349
 350         // look for tag in this alignment and get tag type
 351         char tagType(0);
 352         if ( !al.GetTagType(m_settings->TagToSplit, tagType) )
 353             continue;
 354
 355         // request split method based on tag type
 356         // pass it the current alignment found
 357         switch ( tagType ) {
 358
 359             case (Constants::BAM_TAG_TYPE_INT8)  :
 360             case (Constants::BAM_TAG_TYPE_INT16) :
 361             case (Constants::BAM_TAG_TYPE_INT32) :
 362                 return SplitTagImpl<int32_t>(al);
 363
 364             case (Constants::BAM_TAG_TYPE_UINT8)  :
 365             case (Constants::BAM_TAG_TYPE_UINT16) :
 366             case (Constants::BAM_TAG_TYPE_UINT32) :
 367                 return SplitTagImpl<uint32_t>(al);
 368
 369             case (Constants::BAM_TAG_TYPE_FLOAT)  :
 370                 return SplitTagImpl<float>(al);
 371
 372             case (Constants::BAM_TAG_TYPE_ASCII)  :
 373             case (Constants::BAM_TAG_TYPE_STRING) :
 374             case (Constants::BAM_TAG_TYPE_HEX)    :
 375                 return SplitTagImpl<string>(al);
 376
 377             default:
 378                 fprintf(stderr, "bamtools split ERROR: unknown tag type encountered: [%c]\n", tagType);
 379                 return false;
 380         }
 381     }
 382
 383     // tag not found, but that's not an error - return success
 384     return true;
 385 }
 386
 387 // --------------------------------------------------------------------------------
 388 // template method implementation
 389 // *Technical Note* - use of template methods declared & defined in ".cpp" file
 390 //                    goes against normal practices, but works here because these
 391 //                    are purely internal (no one can call from outside this file)
 392
 393 // close BamWriters & delete pointers
 394 template<typename T>
 395 void SplitTool::SplitToolPrivate::CloseWriters(map<T, BamWriter*>& writers) {
 396
 397     typedef map<T, BamWriter*> WriterMap;
 398     typedef typename WriterMap::iterator WriterMapIterator;
 399
 400     // iterate over writers
 401     WriterMapIterator writerIter = writers.begin();
 402     WriterMapIterator writerEnd  = writers.end();
 403     for ( ; writerIter != writerEnd; ++writerIter ) {
 404         BamWriter* writer = (*writerIter).second;
 405         if ( writer == 0 ) continue;
 406
 407         // close BamWriter
 408         writer->Close();
 409
 410         // destroy BamWriter
 411         delete writer;
 412         writer = 0;
 413     }
 414
 415     // clear the container (destroying the items doesn't remove them)
 416     writers.clear();
 417 }
 418
 419 // handle the various types that are possible for tags
 420 template<typename T>
 421 bool SplitTool::SplitToolPrivate::SplitTagImpl(BamAlignment& al) {
 422
 423     typedef T TagValueType;
 424     typedef map<TagValueType, BamWriter*> WriterMap;
 425     typedef typename WriterMap::iterator WriterMapIterator;
 426
 427     // set up splitting data structure
 428     WriterMap outputFiles;
 429     WriterMapIterator writerIter;
 430
 431     // local variables
 432     const string tag = m_settings->TagToSplit;
 433     BamWriter* writer;
 434     stringstream outputFilenameStream("");
 435     TagValueType currentValue;
 436
 437     // retrieve first alignment tag value
 438     if ( al.GetTag(tag, currentValue) ) {
 439
 440         // open new BamWriter, save first alignment
 441         outputFilenameStream << m_outputFilenameStub << ".TAG_" << tag << "_" << currentValue << ".bam";
 442         writer = new BamWriter;
 443         if ( !writer->Open(outputFilenameStream.str(), m_header, m_references) ) {
 444             cerr << "bamtools split ERROR: could not open " << outputFilenameStream.str()
 445                  << " for writing." << endl;
 446             return false;
 447         }
 448         writer->SaveAlignment(al);
 449
 450         // store in map
 451         outputFiles.insert( make_pair(currentValue, writer) );
 452
 453         // reset stream
 454         outputFilenameStream.str("");
 455     }
 456
 457     // iterate through remaining alignments
 458     while ( m_reader.GetNextAlignment(al) ) {
 459
 460         // skip if this alignment doesn't have TAG
 461         if ( !al.GetTag(tag, currentValue) ) continue;
 462
 463         // look up tag value in map
 464         writerIter = outputFiles.find(currentValue);
 465
 466         // if no writer associated with this value
 467         if ( writerIter == outputFiles.end() ) {
 468
 469             // open new BamWriter
 470             outputFilenameStream << m_outputFilenameStub << ".TAG_" << tag << "_" << currentValue << ".bam";
 471             writer = new BamWriter;
 472             if ( !writer->Open(outputFilenameStream.str(), m_header, m_references) ) {
 473                 cerr << "bamtool split ERROR: could not open " << outputFilenameStream.str()
 474                      << " for writing." << endl;
 475                 return false;
 476             }
 477
 478             // store in map
 479             outputFiles.insert( make_pair(currentValue, writer) );
 480
 481             // reset stream
 482             outputFilenameStream.str("");
 483         }
 484
 485         // else grab corresponding writer
 486         else writer = (*writerIter).second;
 487
 488         // store alignment in proper BAM output file
 489         if ( writer )
 490             writer->SaveAlignment(al);
 491     }
 492
 493     // clean up BamWriters
 494     CloseWriters(outputFiles);
 495
 496     // return success
 497     return true;
 498 }
 499
 500 // ---------------------------------------------
 501 // SplitTool implementation
 502
 503 SplitTool::SplitTool(void)
 504     : AbstractTool()
 505     , m_settings(new SplitSettings)
 506     , m_impl(0)
 507 {
 508     // set program details
 509     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> > ");
 510
 511     // set up options
 512     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 513     Options::AddValueOption("-in",   "BAM filename",  "the input BAM file",  "", m_settings->HasInputFilename,  m_settings->InputFilename,  IO_Opts, Options::StandardIn());
 514     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);
 515
 516     OptionGroup* SplitOpts = Options::CreateOptionGroup("Split Options");
 517     Options::AddOption("-mapped",    "split mapped/unmapped alignments",       m_settings->IsSplittingMapped,    SplitOpts);
 518     Options::AddOption("-paired",    "split single-end/paired-end alignments", m_settings->IsSplittingPaired,    SplitOpts);
 519     Options::AddOption("-reference", "split alignments by reference",          m_settings->IsSplittingReference, SplitOpts);
 520     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)", "",
 521                             m_settings->IsSplittingTag, m_settings->TagToSplit, SplitOpts);
 522 }
 523
 524 SplitTool::~SplitTool(void) {
 525
 526     delete m_settings;
 527     m_settings = 0;
 528
 529     delete m_impl;
 530     m_impl = 0;
 531 }
 532
 533 int SplitTool::Help(void) {
 534     Options::DisplayHelp();
 535     return 0;
 536 }
 537
 538 int SplitTool::Run(int argc, char* argv[]) {
 539
 540     // parse command line arguments
 541     Options::Parse(argc, argv, 1);
 542
 543     // initialize SplitTool with settings
 544     m_impl = new SplitToolPrivate(m_settings);
 545
 546     // run SplitTool, return success/fail
 547     if ( m_impl->Run() )
 548         return 0;
 549     else
 550         return 1;
 551 }