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