src/toolkit/bamtools_count.cpp

   1 // ***************************************************************************
   2 // bamtools_count.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 7 April 2011
   6 // ---------------------------------------------------------------------------
   7 // Prints alignment count for BAM file(s)
   8 // ***************************************************************************
   9
  10 #include "bamtools_count.h"
  11
  12 #include <api/BamMultiReader.h>
  13 #include <utils/bamtools_options.h>
  14 #include <utils/bamtools_utilities.h>
  15 using namespace BamTools;
  16
  17 #include <iostream>
  18 #include <string>
  19 #include <vector>
  20 using namespace std;
  21
  22 // ---------------------------------------------
  23 // CountSettings implementation
  24
  25 struct CountTool::CountSettings {
  26
  27     // flags
  28     bool HasInput;
  29     bool HasRegion;
  30
  31     // filenames
  32     vector<string> InputFiles;
  33     string Region;
  34
  35     // constructor
  36     CountSettings(void)
  37         : HasInput(false)
  38         , HasRegion(false)
  39     { }
  40 };
  41
  42 // ---------------------------------------------
  43 // CountToolPrivate implementation
  44
  45 struct CountTool::CountToolPrivate {
  46
  47     // ctor & dtro
  48     public:
  49         CountToolPrivate(CountTool::CountSettings* settings)
  50             : m_settings(settings)
  51         { }
  52
  53         ~CountToolPrivate(void) { }
  54
  55     // interface
  56     public:
  57         bool Run(void);
  58
  59     // data members
  60     private:
  61         CountTool::CountSettings* m_settings;
  62 };
  63
  64 bool CountTool::CountToolPrivate::Run(void) {
  65
  66     // if no '-in' args supplied, default to stdin
  67     if ( !m_settings->HasInput )
  68         m_settings->InputFiles.push_back(Options::StandardIn());
  69
  70     // open reader without index
  71     BamMultiReader reader;
  72     if ( !reader.Open(m_settings->InputFiles) ) {
  73         cerr << "bamtools count ERROR: could not open input BAM file(s)... Aborting." << endl;
  74         return false;
  75     }
  76
  77     // alignment counter
  78     BamAlignment al;
  79     int alignmentCount(0);
  80
  81     // if no region specified, count entire file
  82     if ( !m_settings->HasRegion ) {
  83         while ( reader.GetNextAlignmentCore(al) )
  84             ++alignmentCount;
  85     }
  86
  87     // otherwise attempt to use region as constraint
  88     else {
  89
  90         // if region string parses OK
  91         BamRegion region;
  92         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
  93
  94             // attempt to find index files
  95             reader.LocateIndexes();
  96
  97             // if index data available for all BAM files, we can use SetRegion
  98             if ( reader.IsIndexLoaded() ) {
  99
 100                 // attempt to set region on reader
 101                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 102                     cerr << "bamtools count ERROR: set region failed. Check that REGION describes a valid range" << endl;
 103                     reader.Close();
 104                     return false;
 105                 }
 106
 107                 // everything checks out, just iterate through specified region, counting alignments
 108                 while ( reader.GetNextAlignmentCore(al) )
 109                     ++alignmentCount;
 110             }
 111
 112             // no index data available, we have to iterate through until we
 113             // find overlapping alignments
 114             else {
 115                 while ( reader.GetNextAlignmentCore(al) ) {
 116                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 117                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 118                     {
 119                         ++alignmentCount;
 120                     }
 121                 }
 122             }
 123         }
 124
 125         // error parsing REGION string
 126         else {
 127             cerr << "bamtools count ERROR: could not parse REGION - " << m_settings->Region << endl;
 128             cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 129                  << endl;
 130             reader.Close();
 131             return false;
 132         }
 133     }
 134
 135     // print results
 136     cout << alignmentCount << endl;
 137
 138     // clean up & exit
 139     reader.Close();
 140     return true;
 141 }
 142
 143 // ---------------------------------------------
 144 // CountTool implementation
 145
 146 CountTool::CountTool(void)
 147     : AbstractTool()
 148     , m_settings(new CountSettings)
 149     , m_impl(0)
 150 {
 151     // set program details
 152     Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)", "[-in <filename> -in <filename> ...] [-region <REGION>]");
 153
 154     // set up options
 155     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 156     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInput,  m_settings->InputFiles, IO_Opts, Options::StandardIn());
 157     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);
 158 }
 159
 160 CountTool::~CountTool(void) {
 161
 162     delete m_settings;
 163     m_settings = 0;
 164
 165     delete m_impl;
 166     m_impl = 0;
 167 }
 168
 169 int CountTool::Help(void) {
 170     Options::DisplayHelp();
 171     return 0;
 172 }
 173
 174 int CountTool::Run(int argc, char* argv[]) {
 175
 176     // parse command line arguments
 177     Options::Parse(argc, argv, 1);
 178
 179     // initialize CountTool with settings
 180     m_impl = new CountToolPrivate(m_settings);
 181
 182     // run CountTool, return success/fail
 183     if ( m_impl->Run() )
 184         return 0;
 185     else
 186         return 1;
 187 }