src/toolkit/bamtools_merge.cpp

   1 // ***************************************************************************
   2 // bamtools_merge.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 7 April 2011
   6 // ---------------------------------------------------------------------------
   7 // Merges multiple BAM files into one
   8 // ***************************************************************************
   9
  10 #include "bamtools_merge.h"
  11
  12 #include <api/BamMultiReader.h>
  13 #include <api/BamWriter.h>
  14 #include <utils/bamtools_options.h>
  15 #include <utils/bamtools_utilities.h>
  16 using namespace BamTools;
  17
  18 #include <iostream>
  19 #include <string>
  20 #include <vector>
  21 using namespace std;
  22
  23 // ---------------------------------------------
  24 // MergeSettings implementation
  25
  26 struct MergeTool::MergeSettings {
  27
  28     // flags
  29     bool HasInputBamFilename;
  30     bool HasOutputBamFilename;
  31     bool IsForceCompression;
  32     bool HasRegion;
  33
  34     // filenames
  35     vector<string> InputFiles;
  36
  37     // other parameters
  38     string OutputFilename;
  39     string Region;
  40
  41     // constructor
  42     MergeSettings(void)
  43         : HasInputBamFilename(false)
  44         , HasOutputBamFilename(false)
  45         , IsForceCompression(false)
  46         , HasRegion(false)
  47         , OutputFilename(Options::StandardOut())
  48     { }
  49 };
  50
  51 // ---------------------------------------------
  52 // MergeToolPrivate implementation
  53
  54 struct MergeTool::MergeToolPrivate {
  55
  56     // ctor & dtor
  57     public:
  58         MergeToolPrivate(MergeTool::MergeSettings* settings)
  59             : m_settings(settings)
  60         { }
  61
  62         ~MergeToolPrivate(void) { }
  63
  64     // interface
  65     public:
  66         bool Run(void);
  67
  68     // data members
  69     private:
  70         MergeTool::MergeSettings* m_settings;
  71 };
  72
  73 bool MergeTool::MergeToolPrivate::Run(void) {
  74
  75     // set to default input if none provided
  76     if ( !m_settings->HasInputBamFilename )
  77         m_settings->InputFiles.push_back(Options::StandardIn());
  78
  79     // opens the BAM files (by default without checking for indexes)
  80     BamMultiReader reader;
  81     if ( !reader.Open(m_settings->InputFiles) ) {
  82         cerr << "bamtools merge ERROR: could not open input BAM file(s)... Aborting." << endl;
  83         return false;
  84     }
  85
  86     // retrieve header & reference dictionary info
  87     std::string mergedHeader = reader.GetHeaderText();
  88     RefVector references = reader.GetReferenceData();
  89
  90     // determine compression mode for BamWriter
  91     bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() &&
  92                                !m_settings->IsForceCompression );
  93     BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
  94     if ( writeUncompressed ) compressionMode = BamWriter::Uncompressed;
  95
  96     // open BamWriter
  97     BamWriter writer;
  98     writer.SetCompressionMode(compressionMode);
  99     if ( !writer.Open(m_settings->OutputFilename, mergedHeader, references) ) {
 100         cerr << "bamtools merge ERROR: could not open "
 101              << m_settings->OutputFilename << " for writing." << endl;
 102         reader.Close();
 103         return false;
 104     }
 105
 106     // if no region specified, store entire contents of file(s)
 107     if ( !m_settings->HasRegion ) {
 108         BamAlignment al;
 109         while ( reader.GetNextAlignmentCore(al) )
 110             writer.SaveAlignment(al);
 111     }
 112
 113     // otherwise attempt to use region as constraint
 114     else {
 115
 116         // if region string parses OK
 117         BamRegion region;
 118         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 119
 120             // attempt to find index files
 121             reader.LocateIndexes();
 122
 123             // if index data available for all BAM files, we can use SetRegion
 124             if ( reader.HasIndexes() ) {
 125
 126                 // attempt to use SetRegion(), if failed report error
 127                 if ( !reader.SetRegion(region.LeftRefID,
 128                                        region.LeftPosition,
 129                                        region.RightRefID,
 130                                        region.RightPosition) )
 131                 {
 132                     cerr << "bamtools merge ERROR: set region failed. Check that REGION describes a valid range"
 133                          << endl;
 134                     reader.Close();
 135                     return false;
 136                 }
 137
 138                 // everything checks out, just iterate through specified region, storing alignments
 139                 BamAlignment al;
 140                 while ( reader.GetNextAlignmentCore(al) )
 141                     writer.SaveAlignment(al);
 142             }
 143
 144             // no index data available, we have to iterate through until we
 145             // find overlapping alignments
 146             else {
 147                 BamAlignment al;
 148                 while ( reader.GetNextAlignmentCore(al) ) {
 149                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 150                          (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 151                     {
 152                         writer.SaveAlignment(al);
 153                     }
 154                 }
 155             }
 156         }
 157
 158         // error parsing REGION string
 159         else {
 160             cerr << "bamtools merge ERROR: could not parse REGION - " << m_settings->Region << endl;
 161             cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 162                  << endl;
 163             reader.Close();
 164             writer.Close();
 165             return false;
 166         }
 167     }
 168
 169     // clean & exit
 170     reader.Close();
 171     writer.Close();
 172     return true;
 173 }
 174
 175 // ---------------------------------------------
 176 // MergeTool implementation
 177
 178 MergeTool::MergeTool(void)
 179     : AbstractTool()
 180     , m_settings(new MergeSettings)
 181     , m_impl(0)
 182 {
 183     // set program details
 184     Options::SetProgramInfo("bamtools merge", "merges multiple BAM files into one", "[-in <filename> -in <filename> ...] [-out <filename> | [-forceCompression]] [-region <REGION>]");
 185
 186     // set up options
 187     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 188     Options::AddValueOption("-in",  "BAM filename", "the input BAM file(s)", "", m_settings->HasInputBamFilename,  m_settings->InputFiles,     IO_Opts);
 189     Options::AddValueOption("-out", "BAM filename", "the output BAM file",   "", m_settings->HasOutputBamFilename, m_settings->OutputFilename, IO_Opts);
 190     Options::AddOption("-forceCompression", "if results are sent to stdout (like when piping to another tool), default behavior is to leave output uncompressed. Use this flag to override and force compression", m_settings->IsForceCompression, IO_Opts);
 191     Options::AddValueOption("-region", "REGION", "genomic region. See README for more details", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
 192 }
 193
 194 MergeTool::~MergeTool(void) {
 195
 196     delete m_settings;
 197     m_settings = 0;
 198
 199     delete m_impl;
 200     m_impl = 0;
 201 }
 202
 203 int MergeTool::Help(void) {
 204     Options::DisplayHelp();
 205     return 0;
 206 }
 207
 208 int MergeTool::Run(int argc, char* argv[]) {
 209
 210     // parse command line arguments
 211     Options::Parse(argc, argv, 1);
 212
 213     // initialize MergeTool with settings
 214     m_impl = new MergeToolPrivate(m_settings);
 215
 216     // run MergeTool, return success/fail
 217     if ( m_impl->Run() )
 218         return 0;
 219     else
 220         return 1;
 221 }