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