1 // ***************************************************************************
2 // bamtools_count.cpp (c) 2010 Derek Barnett, Erik Garrison
3 // Marth Lab, Department of Biology, Boston College
4 // ---------------------------------------------------------------------------
5 // Last modified: 10 December 2012
6 // ---------------------------------------------------------------------------
7 // Prints alignment count for BAM file(s)
8 // ***************************************************************************
10 #include "bamtools_count.h"
12 #include <api/BamAlgorithms.h>
13 #include <api/BamMultiReader.h>
14 #include <utils/bamtools_options.h>
15 #include <utils/bamtools_utilities.h>
16 using namespace BamTools;
24 // ---------------------------------------------
25 // CountSettings implementation
27 struct CountTool::CountSettings {
31 bool HasInputFilelist;
35 vector<string> InputFiles;
42 , HasInputFilelist(false)
47 // ---------------------------------------------
48 // CountToolPrivate implementation
50 struct CountTool::CountToolPrivate {
54 CountToolPrivate(CountTool::CountSettings* settings)
55 : m_settings(settings)
58 ~CountToolPrivate(void) { }
66 CountTool::CountSettings* m_settings;
69 bool CountTool::CountToolPrivate::Run(void) {
71 // set to default input if none provided
72 if ( !m_settings->HasInput && !m_settings->HasInputFilelist )
73 m_settings->InputFiles.push_back(Options::StandardIn());
75 // add files in the filelist to the input file list
76 if ( m_settings->HasInputFilelist ) {
78 ifstream filelist(m_settings->InputFilelist.c_str(), ios::in);
79 if ( !filelist.is_open() ) {
80 cerr << "bamtools count ERROR: could not open input BAM file list... Aborting." << endl;
85 while ( getline(filelist, line) )
86 m_settings->InputFiles.push_back(line);
89 // open reader without index
90 BamMultiReader reader;
91 if ( !reader.Open(m_settings->InputFiles) ) {
92 cerr << "bamtools count ERROR: could not open input BAM file(s)... Aborting." << endl;
98 int alignmentCount(0);
100 // if no region specified, count entire file
101 if ( !m_settings->HasRegion ) {
102 while ( reader.GetNextAlignmentCore(al) )
106 // otherwise attempt to use region as constraint
109 // if region string parses OK
111 if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
113 // attempt to find index files
114 reader.LocateIndexes();
116 // if index data available for all BAM files, we can use SetRegion
117 if ( reader.HasIndexes() ) {
119 // attempt to set region on reader
120 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
121 cerr << "bamtools count ERROR: set region failed. Check that REGION describes a valid range" << endl;
126 // everything checks out, just iterate through specified region, counting alignments
127 while ( reader.GetNextAlignmentCore(al) )
131 // no index data available, we have to iterate through until we
132 // find overlapping alignments
134 while ( reader.GetNextAlignmentCore(al) ) {
135 if ( (al.RefID >= region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) &&
136 (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
144 // error parsing REGION string
146 cerr << "bamtools count ERROR: could not parse REGION - " << m_settings->Region << endl;
147 cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
155 cout << alignmentCount << endl;
162 // ---------------------------------------------
163 // CountTool implementation
165 CountTool::CountTool(void)
167 , m_settings(new CountSettings)
170 // set program details
171 Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)",
172 "[-in <filename> -in <filename> ... | -list <filelist>] [-region <REGION>]");
175 OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
176 Options::AddValueOption("-in", "BAM filename", "the input BAM file(s)", "", m_settings->HasInput, m_settings->InputFiles, IO_Opts, Options::StandardIn());
177 Options::AddValueOption("-list", "filename", "the input BAM file list, one line per file", "", m_settings->HasInputFilelist, m_settings->InputFilelist, IO_Opts);
178 Options::AddValueOption("-region", "REGION",
179 "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",
180 "", m_settings->HasRegion, m_settings->Region, IO_Opts);
183 CountTool::~CountTool(void) {
192 int CountTool::Help(void) {
193 Options::DisplayHelp();
197 int CountTool::Run(int argc, char* argv[]) {
199 // parse command line arguments
200 Options::Parse(argc, argv, 1);
202 // initialize CountTool with settings
203 m_impl = new CountToolPrivate(m_settings);
205 // run CountTool, return success/fail