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