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/BamAlgorithms.h>
  13 #include <api/BamMultiReader.h>
  14 #include <utils/bamtools_options.h>
  15 #include <utils/bamtools_utilities.h>
  16 using namespace BamTools;
  17
  18 #include <iostream>
  19 #include <string>
  20 #include <vector>
  21 using namespace std;
  22
  23 // ---------------------------------------------
  24 // CountSettings implementation
  25
  26 struct CountTool::CountSettings {
  27
  28     // flags
  29     bool HasInput;
  30     bool HasRegion;
  31
  32     // filenames
  33     vector<string> InputFiles;
  34     string Region;
  35
  36     // constructor
  37     CountSettings(void)
  38         : HasInput(false)
  39         , HasRegion(false)
  40     { }
  41 };
  42
  43 // ---------------------------------------------
  44 // CountToolPrivate implementation
  45
  46 struct CountTool::CountToolPrivate {
  47
  48     // ctor & dtro
  49     public:
  50         CountToolPrivate(CountTool::CountSettings* settings)
  51             : m_settings(settings)
  52         { }
  53
  54         ~CountToolPrivate(void) { }
  55
  56     // interface
  57     public:
  58         bool Run(void);
  59
  60     // data members
  61     private:
  62         CountTool::CountSettings* m_settings;
  63 };
  64
  65 bool CountTool::CountToolPrivate::Run(void) {
  66
  67     // if no '-in' args supplied, default to stdin
  68     if ( !m_settings->HasInput )
  69         m_settings->InputFiles.push_back(Options::StandardIn());
  70
  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;
  75         return false;
  76     }
  77
  78     // alignment counter
  79     BamAlignment al;
  80     int alignmentCount(0);
  81
  82     // if no region specified, count entire file
  83     if ( !m_settings->HasRegion ) {
  84         while ( reader.GetNextAlignmentCore(al) )
  85             ++alignmentCount;
  86     }
  87
  88     // otherwise attempt to use region as constraint
  89     else {
  90
  91         // if region string parses OK
  92         BamRegion region;
  93         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
  94
  95             // attempt to find index files
  96             reader.LocateIndexes();
  97
  98             // if index data available for all BAM files, we can use SetRegion
  99             if ( reader.HasIndexes() ) {
 100
 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;
 104                     reader.Close();
 105                     return false;
 106                 }
 107
 108                 // everything checks out, just iterate through specified region, counting alignments
 109                 while ( reader.GetNextAlignmentCore(al) )
 110                     ++alignmentCount;
 111             }
 112
 113             // no index data available, we have to iterate through until we
 114             // find overlapping alignments
 115             else {
 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) )
 119                     {
 120                         ++alignmentCount;
 121                     }
 122                 }
 123             }
 124         }
 125
 126         // error parsing REGION string
 127         else {
 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"
 130                  << endl;
 131             reader.Close();
 132             return false;
 133         }
 134     }
 135
 136     // print results
 137     cout << alignmentCount << endl;
 138
 139     // clean up & exit
 140     reader.Close();
 141     return true;
 142 }
 143
 144 // ---------------------------------------------
 145 // CountTool implementation
 146
 147 CountTool::CountTool(void)
 148     : AbstractTool()
 149     , m_settings(new CountSettings)
 150     , m_impl(0)
 151 {
 152     // set program details
 153     Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)", "[-in <filename> -in <filename> ...] [-region <REGION>]");
 154
 155     // set up options
 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);
 159 }
 160
 161 CountTool::~CountTool(void) {
 162
 163     delete m_settings;
 164     m_settings = 0;
 165
 166     delete m_impl;
 167     m_impl = 0;
 168 }
 169
 170 int CountTool::Help(void) {
 171     Options::DisplayHelp();
 172     return 0;
 173 }
 174
 175 int CountTool::Run(int argc, char* argv[]) {
 176
 177     // parse command line arguments
 178     Options::Parse(argc, argv, 1);
 179
 180     // initialize CountTool with settings
 181     m_impl = new CountToolPrivate(m_settings);
 182
 183     // run CountTool, return success/fail
 184     if ( m_impl->Run() )
 185         return 0;
 186     else
 187         return 1;
 188 }