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 September 2010
   7 // ---------------------------------------------------------------------------
   8 // Merges multiple BAM files into one.
   9 // ***************************************************************************
  10
  11 #include <iostream>
  12 #include <string>
  13 #include <vector>
  14
  15 #include "bamtools_merge.h"
  16 #include "bamtools_options.h"
  17 #include "bamtools_utilities.h"
  18 #include "BamMultiReader.h"
  19 #include "BamWriter.h"
  20
  21 using namespace std;
  22 using namespace BamTools;
  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 // MergeTool implementation
  54
  55 MergeTool::MergeTool(void)
  56     : AbstractTool()
  57     , m_settings(new MergeSettings)
  58 {
  59     // set program details
  60     Options::SetProgramInfo("bamtools merge", "merges multiple BAM files into one", "[-in <filename> -in <filename> ...] [-out <filename>]");
  61
  62     // set up options
  63     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  64     Options::AddValueOption("-in",  "BAM filename", "the input BAM file(s)", "", m_settings->HasInputBamFilename,  m_settings->InputFiles,     IO_Opts);
  65     Options::AddValueOption("-out", "BAM filename", "the output BAM file",   "", m_settings->HasOutputBamFilename, m_settings->OutputFilename, IO_Opts);
  66     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);
  67
  68     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  69     Options::AddValueOption("-region", "REGION", "genomic region. See README for more details", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
  70 }
  71
  72 MergeTool::~MergeTool(void) {
  73     delete m_settings;
  74     m_settings = 0;
  75 }
  76
  77 int MergeTool::Help(void) {
  78     Options::DisplayHelp();
  79     return 0;
  80 }
  81
  82 int MergeTool::Run(int argc, char* argv[]) {
  83
  84     // parse command line arguments
  85     Options::Parse(argc, argv, 1);
  86
  87      // set to default input if none provided
  88     if ( !m_settings->HasInputBamFilename )
  89         m_settings->InputFiles.push_back(Options::StandardIn());
  90
  91     // opens the BAM files (by default without checking for indexes)
  92     BamMultiReader reader;
  93     if ( !reader.Open(m_settings->InputFiles, false, true) ) {
  94         cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
  95         return 1;
  96     }
  97
  98     // retrieve header & reference dictionary info
  99     std::string mergedHeader = reader.GetHeaderText();
 100     RefVector references = reader.GetReferenceData();
 101
 102     // open writer
 103     BamWriter writer;
 104     bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() && !m_settings->IsForceCompression );
 105     if ( !writer.Open(m_settings->OutputFilename, mergedHeader, references, writeUncompressed) ) {
 106         cerr << "ERROR: Could not open BAM file " << m_settings->OutputFilename << " for writing... Aborting." << endl;
 107         reader.Close();
 108         return 1;
 109     }
 110
 111     // if no region specified, store entire contents of file(s)
 112     if ( !m_settings->HasRegion ) {
 113         BamAlignment al;
 114         while ( reader.GetNextAlignmentCore(al) )
 115             writer.SaveAlignment(al);
 116     }
 117
 118     // otherwise attempt to use region as constraint
 119     else {
 120
 121         // if region string parses OK
 122         BamRegion region;
 123         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 124
 125             // attempt to re-open reader with index files
 126             reader.Close();
 127             bool openedOK = reader.Open(m_settings->InputFiles, true, true );
 128
 129             // if error
 130             if ( !openedOK ) {
 131                 cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
 132                 return 1;
 133             }
 134
 135             // if index data available, we can use SetRegion
 136             if ( reader.IsIndexLoaded() ) {
 137
 138                 // attempt to use SetRegion(), if failed report error
 139                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 140                     cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
 141                     reader.Close();
 142                     return 1;
 143                 }
 144
 145                 // everything checks out, just iterate through specified region, storing alignments
 146                 BamAlignment al;
 147                 while ( reader.GetNextAlignmentCore(al) )
 148                     writer.SaveAlignment(al);
 149             }
 150
 151             // no index data available, we have to iterate through until we
 152             // find overlapping alignments
 153             else {
 154                 BamAlignment al;
 155                 while ( reader.GetNextAlignmentCore(al) ) {
 156                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 157                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 158                     {
 159                         writer.SaveAlignment(al);
 160                     }
 161                 }
 162             }
 163         }
 164
 165         // error parsing REGION string
 166         else {
 167             cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
 168             cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 169             reader.Close();
 170             writer.Close();
 171             return 1;
 172         }
 173     }
 174
 175     // clean & exit
 176     reader.Close();
 177     writer.Close();
 178     return 0;
 179 }