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 // ***************************************************************************
11 #include "bamtools_coverage.h"
13 #include <api/BamReader.h>
14 #include <utils/bamtools_options.h>
15 #include <utils/bamtools_pileup_engine.h>
16 using namespace BamTools;
26 // ---------------------------------------------
27 // CoverageVisitor implementation
29 class CoverageVisitor : public PileupVisitor {
32 CoverageVisitor(const RefVector& references, ostream* out)
34 , m_references(references)
37 ~CoverageVisitor(void) { }
39 // PileupVisitor interface implementation
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;
49 RefVector m_references;
53 } // namespace BamTools
55 // ---------------------------------------------
56 // CoverageSettings implementation
58 struct CoverageTool::CoverageSettings {
65 string InputBamFilename;
66 string OutputFilename;
69 CoverageSettings(void)
71 , HasOutputFile(false)
72 , InputBamFilename(Options::StandardIn())
73 , OutputFilename(Options::StandardOut())
77 // ---------------------------------------------
78 // CoverageToolPrivate implementation
80 struct CoverageTool::CoverageToolPrivate {
84 CoverageToolPrivate(CoverageTool::CoverageSettings* settings);
85 ~CoverageToolPrivate(void);
93 CoverageTool::CoverageSettings* m_settings;
95 RefVector m_references;
98 CoverageTool::CoverageToolPrivate::CoverageToolPrivate(CoverageTool::CoverageSettings* settings)
99 : m_settings(settings)
100 , m_out(cout.rdbuf()) // default output to cout
103 CoverageTool::CoverageToolPrivate::~CoverageToolPrivate(void) { }
105 bool CoverageTool::CoverageToolPrivate::Run(void) {
107 // if output filename given
109 if ( m_settings->HasOutputFile ) {
111 // open output file stream
112 outFile.open(m_settings->OutputFilename.c_str());
114 cerr << "bamtools coverage ERROR: could not open " << m_settings->OutputFilename
115 << " for output" << endl;
119 // set m_out to file's streambuf
120 m_out.rdbuf(outFile.rdbuf());
123 //open our BAM reader
125 if ( !reader.Open(m_settings->InputBamFilename) ) {
126 cerr << "bamtools coverage ERROR: could not open input BAM file: " << m_settings->InputBamFilename << endl;
130 // retrieve references
131 m_references = reader.GetReferenceData();
133 // set up our output 'visitor'
134 CoverageVisitor* cv = new CoverageVisitor(m_references, &m_out);
136 // set up pileup engine with 'visitor'
138 pileup.AddVisitor(cv);
140 // process input data
142 while ( reader.GetNextAlignment(al) )
143 pileup.AddAlignment(al);
147 if ( m_settings->HasOutputFile )
156 // ---------------------------------------------
157 // CoverageTool implementation
159 CoverageTool::CoverageTool(void)
161 , m_settings(new CoverageSettings)
164 // set program details
165 Options::SetProgramInfo("bamtools coverage", "prints coverage data for a single BAM file", "[-in <filename>] [-out <filename>]");
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());
173 CoverageTool::~CoverageTool(void) {
181 int CoverageTool::Help(void) {
182 Options::DisplayHelp();
186 int CoverageTool::Run(int argc, char* argv[]) {
188 // parse command line arguments
189 Options::Parse(argc, argv, 1);
191 // run internal ConvertTool implementation, return success/fail
192 m_impl = new CoverageToolPrivate(m_settings);