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