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