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: 10 December 2012
   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 <fstream>
  19 #include <iostream>
  20 #include <string>
  21 #include <vector>
  22 using namespace std;
  23
  24 // ---------------------------------------------
  25 // CountSettings implementation
  26
  27 struct CountTool::CountSettings {
  28
  29     // flags
  30     bool HasInput;
  31     bool HasInputFilelist;
  32     bool HasRegion;
  33
  34     // filenames
  35     vector<string> InputFiles;
  36     string InputFilelist;
  37     string Region;
  38
  39     // constructor
  40     CountSettings(void)
  41         : HasInput(false)
  42         , HasInputFilelist(false)
  43         , HasRegion(false)
  44     { }
  45 };
  46
  47 // ---------------------------------------------
  48 // CountToolPrivate implementation
  49
  50 struct CountTool::CountToolPrivate {
  51
  52     // ctor & dtro
  53     public:
  54         CountToolPrivate(CountTool::CountSettings* settings)
  55             : m_settings(settings)
  56         { }
  57
  58         ~CountToolPrivate(void) { }
  59
  60     // interface
  61     public:
  62         bool Run(void);
  63
  64     // data members
  65     private:
  66         CountTool::CountSettings* m_settings;
  67 };
  68
  69 bool CountTool::CountToolPrivate::Run(void) {
  70
  71     // set to default input if none provided
  72     if ( !m_settings->HasInput && !m_settings->HasInputFilelist )
  73         m_settings->InputFiles.push_back(Options::StandardIn());
  74
  75     // add files in the filelist to the input file list
  76     if ( m_settings->HasInputFilelist ) {
  77
  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;
  81             return false;
  82         }
  83
  84         string line;
  85         while ( getline(filelist, line) )
  86             m_settings->InputFiles.push_back(line);
  87     }
  88
  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;
  93         return false;
  94     }
  95
  96     // alignment counter
  97     BamAlignment al;
  98     int alignmentCount(0);
  99
 100     // if no region specified, count entire file
 101     if ( !m_settings->HasRegion ) {
 102         while ( reader.GetNextAlignmentCore(al) )
 103             ++alignmentCount;
 104     }
 105
 106     // otherwise attempt to use region as constraint
 107     else {
 108
 109         // if region string parses OK
 110         BamRegion region;
 111         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 112
 113             // attempt to find index files
 114             reader.LocateIndexes();
 115
 116             // if index data available for all BAM files, we can use SetRegion
 117             if ( reader.HasIndexes() ) {
 118
 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;
 122                     reader.Close();
 123                     return false;
 124                 }
 125
 126                 // everything checks out, just iterate through specified region, counting alignments
 127                 while ( reader.GetNextAlignmentCore(al) )
 128                     ++alignmentCount;
 129             }
 130
 131             // no index data available, we have to iterate through until we
 132             // find overlapping alignments
 133             else {
 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) )
 137                     {
 138                         ++alignmentCount;
 139                     }
 140                 }
 141             }
 142         }
 143
 144         // error parsing REGION string
 145         else {
 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"
 148                  << endl;
 149             reader.Close();
 150             return false;
 151         }
 152     }
 153
 154     // print results
 155     cout << alignmentCount << endl;
 156
 157     // clean up & exit
 158     reader.Close();
 159     return true;
 160 }
 161
 162 // ---------------------------------------------
 163 // CountTool implementation
 164
 165 CountTool::CountTool(void)
 166     : AbstractTool()
 167     , m_settings(new CountSettings)
 168     , m_impl(0)
 169 {
 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>]");
 173
 174     // set up options
 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);
 181 }
 182
 183 CountTool::~CountTool(void) {
 184
 185     delete m_settings;
 186     m_settings = 0;
 187
 188     delete m_impl;
 189     m_impl = 0;
 190 }
 191
 192 int CountTool::Help(void) {
 193     Options::DisplayHelp();
 194     return 0;
 195 }
 196
 197 int CountTool::Run(int argc, char* argv[]) {
 198
 199     // parse command line arguments
 200     Options::Parse(argc, argv, 1);
 201
 202     // initialize CountTool with settings
 203     m_impl = new CountToolPrivate(m_settings);
 204
 205     // run CountTool, return success/fail
 206     if ( m_impl->Run() )
 207         return 0;
 208     else
 209         return 1;
 210 }