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