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 //         void Visit(const int& refId, const int& position, const vector<BamAlignment>& alignments) {
  40         void Visit(const PileupPosition& pileupData) {
  41             // -----------------------------------------
  42             // print coverage results ( tab-delimited )
  43             // <refName> <0-based pos> <number of alleles>
  44 //             *m_out << m_references[refId].RefName << "\t" << position << "\t" << alignments.size() << endl;
  45             *m_out << m_references[pileupData.RefId].RefName << "\t"
  46                    << pileupData.Position << "\t"
  47                    << pileupData.PileupAlignments.size() << endl;
  48         }
  49
  50     private:
  51         RefVector m_references;
  52         ostream*  m_out;
  53 };
  54
  55 } // namespace BamTools
  56
  57 // ---------------------------------------------
  58 // CoverageSettings implementation
  59
  60 struct CoverageTool::CoverageSettings {
  61
  62     // flags
  63     bool HasInputFile;
  64     bool HasOutputFile;
  65
  66     // filenames
  67     string InputBamFilename;
  68     string OutputFilename;
  69
  70     // constructor
  71     CoverageSettings(void)
  72         : HasInputFile(false)
  73         , HasOutputFile(false)
  74         , InputBamFilename(Options::StandardIn())
  75         , OutputFilename(Options::StandardOut())
  76     { }
  77 };
  78
  79 // ---------------------------------------------
  80 // CoverageToolPrivate implementation
  81
  82 struct CoverageTool::CoverageToolPrivate {
  83
  84     // ctor & dtor
  85     public:
  86         CoverageToolPrivate(CoverageTool::CoverageSettings* settings);
  87         ~CoverageToolPrivate(void);
  88
  89     // interface
  90     public:
  91         bool Run(void);
  92
  93     // data members
  94     private:
  95         CoverageTool::CoverageSettings* m_settings;
  96         ostream m_out;
  97         RefVector m_references;
  98 };
  99
 100 CoverageTool::CoverageToolPrivate::CoverageToolPrivate(CoverageTool::CoverageSettings* settings)
 101     : m_settings(settings)
 102     , m_out(cout.rdbuf()) // default output to cout
 103 { }
 104
 105 CoverageTool::CoverageToolPrivate::~CoverageToolPrivate(void) { }
 106
 107 bool CoverageTool::CoverageToolPrivate::Run(void) {
 108
 109     // if output filename given
 110     ofstream outFile;
 111     if ( m_settings->HasOutputFile ) {
 112
 113         // open output file stream
 114         outFile.open(m_settings->OutputFilename.c_str());
 115         if ( !outFile ) {
 116             cerr << "Could not open " << m_settings->OutputFilename << " for output." << endl;
 117             return false;
 118         }
 119
 120         // set m_out to file's streambuf
 121         m_out.rdbuf(outFile.rdbuf());
 122     }
 123
 124     //open our BAM reader
 125     BamReader reader;
 126     reader.Open(m_settings->InputBamFilename);
 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
 141     // clean up
 142     reader.Close();
 143     if ( m_settings->HasOutputFile ) outFile.close();
 144     delete cv;
 145     cv = 0;
 146
 147     // return success
 148     return true;
 149 }
 150
 151 // ---------------------------------------------
 152 // CoverageTool implementation
 153
 154 CoverageTool::CoverageTool(void)
 155     : AbstractTool()
 156     , m_settings(new CoverageSettings)
 157     , m_impl(0)
 158 {
 159     // set program details
 160     Options::SetProgramInfo("bamtools coverage", "prints coverage data for a single BAM file", "[-in <filename>] [-out <filename>]");
 161
 162     // set up options
 163     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 164     Options::AddValueOption("-in",  "BAM filename", "the input BAM file", "", m_settings->HasInputFile,  m_settings->InputBamFilename, IO_Opts, Options::StandardIn());
 165     Options::AddValueOption("-out", "filename",     "the output file",    "", m_settings->HasOutputFile, m_settings->OutputFilename,   IO_Opts, Options::StandardOut());
 166 }
 167
 168 CoverageTool::~CoverageTool(void) {
 169     delete m_settings;
 170     m_settings = 0;
 171
 172     delete m_impl;
 173     m_impl = 0;
 174 }
 175
 176 int CoverageTool::Help(void) {
 177     Options::DisplayHelp();
 178     return 0;
 179 }
 180
 181 int CoverageTool::Run(int argc, char* argv[]) {
 182
 183     // parse command line arguments
 184     Options::Parse(argc, argv, 1);
 185
 186     // run internal ConvertTool implementation, return success/fail
 187     m_impl = new CoverageToolPrivate(m_settings);
 188
 189     if ( m_impl->Run() )
 190         return 0;
 191     else
 192         return 1;
 193 }