1 // ***************************************************************************
2 // bamtools_count.cpp (c) 2010 Derek Barnett, Erik Garrison
3 // Marth Lab, Department of Biology, Boston College
4 // All rights reserved.
5 // ---------------------------------------------------------------------------
6 // Last modified: 7 April 2011
7 // ---------------------------------------------------------------------------
8 // Prints alignment count for BAM file(s)
9 // ***************************************************************************
11 #include "bamtools_count.h"
13 #include <api/BamMultiReader.h>
14 #include <utils/bamtools_options.h>
15 #include <utils/bamtools_utilities.h>
16 using namespace BamTools;
23 // ---------------------------------------------
24 // CountSettings implementation
26 struct CountTool::CountSettings {
33 vector<string> InputFiles;
43 // ---------------------------------------------
44 // CountToolPrivate implementation
46 struct CountTool::CountToolPrivate {
50 CountToolPrivate(CountTool::CountSettings* settings)
51 : m_settings(settings)
54 ~CountToolPrivate(void) { }
62 CountTool::CountSettings* m_settings;
65 bool CountTool::CountToolPrivate::Run(void) {
67 // if no '-in' args supplied, default to stdin
68 if ( !m_settings->HasInput )
69 m_settings->InputFiles.push_back(Options::StandardIn());
71 // open reader without index
72 BamMultiReader reader;
73 if ( !reader.Open(m_settings->InputFiles) ) {
74 cerr << "bamtools count ERROR: could not open input BAM file(s)... Aborting." << endl;
80 int alignmentCount(0);
82 // if no region specified, count entire file
83 if ( !m_settings->HasRegion ) {
84 while ( reader.GetNextAlignmentCore(al) )
88 // otherwise attempt to use region as constraint
91 // if region string parses OK
93 if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
95 // attempt to find index files
96 reader.LocateIndexes();
98 // if index data available for all BAM files, we can use SetRegion
99 if ( reader.IsIndexLoaded() ) {
101 // attempt to set region on reader
102 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
103 cerr << "bamtools count ERROR: set region failed. Check that REGION describes a valid range" << endl;
108 // everything checks out, just iterate through specified region, counting alignments
109 while ( reader.GetNextAlignmentCore(al) )
113 // no index data available, we have to iterate through until we
114 // find overlapping alignments
116 while ( reader.GetNextAlignmentCore(al) ) {
117 if ( (al.RefID >= region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) &&
118 (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
126 // error parsing REGION string
128 cerr << "bamtools count ERROR: could not parse REGION - " << m_settings->Region << endl;
129 cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
137 cout << alignmentCount << endl;
144 // ---------------------------------------------
145 // CountTool implementation
147 CountTool::CountTool(void)
149 , m_settings(new CountSettings)
152 // set program details
153 Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)", "[-in <filename> -in <filename> ...] [-region <REGION>]");
156 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
157 Options::AddValueOption("-in", "BAM filename", "the input BAM file(s)", "", m_settings->HasInput, m_settings->InputFiles, IO_Opts, Options::StandardIn());
158 Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance, and is used automatically if it exists. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
161 CountTool::~CountTool(void) {
170 int CountTool::Help(void) {
171 Options::DisplayHelp();
175 int CountTool::Run(int argc, char* argv[]) {
177 // parse command line arguments
178 Options::Parse(argc, argv, 1);
180 // initialize CountTool with settings
181 m_impl = new CountToolPrivate(m_settings);
183 // run CountTool, return success/fail