src/toolkit/bamtools_coverage.cpp

   1 // ***************************************************************************
   2 // bamtools_coverage.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 21 March 2011
   7 // ---------------------------------------------------------------------------
   8 // Prints coverage data for a single BAM file
   9 // ***************************************************************************
  10
  11 #include "bamtools_coverage.h"
  12
  13 #include <api/BamReader.h>
  14 #include <utils/bamtools_options.h>
  15 #include <utils/bamtools_pileup_engine.h>
  16 using namespace BamTools;
  17
  18 #include <iostream>
  19 #include <fstream>
  20 #include <string>
  21 #include <vector>
  22 using namespace std;
  23
  24 namespace BamTools {
  25
  26 // ---------------------------------------------
  27 // CoverageVisitor implementation
  28
  29 class CoverageVisitor : public PileupVisitor {
  30
  31     public:
  32         CoverageVisitor(const RefVector& references, ostream* out)
  33             : PileupVisitor()
  34             , m_references(references)
  35             , m_out(out)
  36         { }
  37         ~CoverageVisitor(void) { }
  38
  39     // PileupVisitor interface implementation
  40     public:
  41         // prints coverage results ( tab-delimited )
  42         void Visit(const PileupPosition& pileupData) {
  43             *m_out << m_references[pileupData.RefId].RefName << "\t"
  44                    << pileupData.Position << "\t"
  45                    << pileupData.PileupAlignments.size() << endl;
  46         }
  47
  48     private:
  49         RefVector m_references;
  50         ostream*  m_out;
  51 };
  52
  53 } // namespace BamTools
  54
  55 // ---------------------------------------------
  56 // CoverageSettings implementation
  57
  58 struct CoverageTool::CoverageSettings {
  59
  60     // flags
  61     bool HasInputFile;
  62     bool HasOutputFile;
  63
  64     // filenames
  65     string InputBamFilename;
  66     string OutputFilename;
  67
  68     // constructor
  69     CoverageSettings(void)
  70         : HasInputFile(false)
  71         , HasOutputFile(false)
  72         , InputBamFilename(Options::StandardIn())
  73         , OutputFilename(Options::StandardOut())
  74     { }
  75 };
  76
  77 // ---------------------------------------------
  78 // CoverageToolPrivate implementation
  79
  80 struct CoverageTool::CoverageToolPrivate {
  81
  82     // ctor & dtor
  83     public:
  84         CoverageToolPrivate(CoverageTool::CoverageSettings* settings);
  85         ~CoverageToolPrivate(void);
  86
  87     // interface
  88     public:
  89         bool Run(void);
  90
  91     // data members
  92     private:
  93         CoverageTool::CoverageSettings* m_settings;
  94         ostream m_out;
  95         RefVector m_references;
  96 };
  97
  98 CoverageTool::CoverageToolPrivate::CoverageToolPrivate(CoverageTool::CoverageSettings* settings)
  99     : m_settings(settings)
 100     , m_out(cout.rdbuf()) // default output to cout
 101 { }
 102
 103 CoverageTool::CoverageToolPrivate::~CoverageToolPrivate(void) { }
 104
 105 bool CoverageTool::CoverageToolPrivate::Run(void) {
 106
 107     // if output filename given
 108     ofstream outFile;
 109     if ( m_settings->HasOutputFile ) {
 110
 111         // open output file stream
 112         outFile.open(m_settings->OutputFilename.c_str());
 113         if ( !outFile ) {
 114             cerr << "bamtools coverage ERROR: could not open " << m_settings->OutputFilename
 115                  << " for output" << endl;
 116             return false;
 117         }
 118
 119         // set m_out to file's streambuf
 120         m_out.rdbuf(outFile.rdbuf());
 121     }
 122
 123     //open our BAM reader
 124     BamReader reader;
 125     if ( !reader.Open(m_settings->InputBamFilename) ) {
 126         cerr << "bamtools coverage ERROR: could not open input BAM file: " << m_settings->InputBamFilename << endl;
 127         return false;
 128     }
 129
 130     // retrieve references
 131     m_references = reader.GetReferenceData();
 132
 133     // set up our output 'visitor'
 134     CoverageVisitor* cv = new CoverageVisitor(m_references, &m_out);
 135
 136     // set up pileup engine with 'visitor'
 137     PileupEngine pileup;
 138     pileup.AddVisitor(cv);
 139
 140     // process input data
 141     BamAlignment al;
 142     while ( reader.GetNextAlignment(al) )
 143         pileup.AddAlignment(al);
 144
 145     // clean up
 146     reader.Close();
 147     if ( m_settings->HasOutputFile )
 148         outFile.close();
 149     delete cv;
 150     cv = 0;
 151
 152     // return success
 153     return true;
 154 }
 155
 156 // ---------------------------------------------
 157 // CoverageTool implementation
 158
 159 CoverageTool::CoverageTool(void)
 160     : AbstractTool()
 161     , m_settings(new CoverageSettings)
 162     , m_impl(0)
 163 {
 164     // set program details
 165     Options::SetProgramInfo("bamtools coverage", "prints coverage data for a single BAM file", "[-in <filename>] [-out <filename>]");
 166
 167     // set up options
 168     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 169     Options::AddValueOption("-in",  "BAM filename", "the input BAM file", "", m_settings->HasInputFile,  m_settings->InputBamFilename, IO_Opts, Options::StandardIn());
 170     Options::AddValueOption("-out", "filename",     "the output file",    "", m_settings->HasOutputFile, m_settings->OutputFilename,   IO_Opts, Options::StandardOut());
 171 }
 172
 173 CoverageTool::~CoverageTool(void) {
 174     delete m_settings;
 175     m_settings = 0;
 176
 177     delete m_impl;
 178     m_impl = 0;
 179 }
 180
 181 int CoverageTool::Help(void) {
 182     Options::DisplayHelp();
 183     return 0;
 184 }
 185
 186 int CoverageTool::Run(int argc, char* argv[]) {
 187
 188     // parse command line arguments
 189     Options::Parse(argc, argv, 1);
 190
 191     // run internal ConvertTool implementation, return success/fail
 192     m_impl = new CoverageToolPrivate(m_settings);
 193
 194     if ( m_impl->Run() )
 195         return 0;
 196     else
 197         return 1;
 198 }