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