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 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 23 March 2011
   7 // ---------------------------------------------------------------------------
   8 // Prints alignment count for BAM file(s)
   9 // ***************************************************************************
  10
  11 #include "bamtools_count.h"
  12
  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 // CountTool implementation
  45
  46 CountTool::CountTool(void)
  47     : AbstractTool()
  48     , m_settings(new CountSettings)
  49 {
  50     // set program details
  51     Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)", "[-in <filename> -in <filename> ...] [-region <REGION>]");
  52
  53     // set up options
  54     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  55     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInput,  m_settings->InputFiles, IO_Opts, Options::StandardIn());
  56     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);
  57 }
  58
  59 CountTool::~CountTool(void) {
  60     delete m_settings;
  61     m_settings = 0;
  62 }
  63
  64 int CountTool::Help(void) {
  65     Options::DisplayHelp();
  66     return 0;
  67 }
  68
  69 int CountTool::Run(int argc, char* argv[]) {
  70
  71     // parse command line arguments
  72     Options::Parse(argc, argv, 1);
  73
  74     // if no '-in' args supplied, default to stdin
  75     if ( !m_settings->HasInput )
  76         m_settings->InputFiles.push_back(Options::StandardIn());
  77
  78     // open reader without index
  79     BamMultiReader reader;
  80     if ( !reader.Open(m_settings->InputFiles) ) {
  81         cerr << "bamtools count ERROR: could not open input BAM file(s)... Aborting." << endl;
  82         return 1;
  83     }
  84
  85     // alignment counter
  86     BamAlignment al;
  87     int alignmentCount(0);
  88
  89     // if no region specified, count entire file
  90     if ( !m_settings->HasRegion ) {
  91         while ( reader.GetNextAlignmentCore(al) )
  92             ++alignmentCount;
  93     }
  94
  95     // otherwise attempt to use region as constraint
  96     else {
  97
  98         // if region string parses OK
  99         BamRegion region;
 100         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 101
 102             // attempt to find index files
 103             reader.LocateIndexes();
 104
 105             // if index data available for all BAM files, we can use SetRegion
 106             if ( reader.IsIndexLoaded() ) {
 107
 108                 // attempt to set region on reader
 109                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 110                     cerr << "bamtools count ERROR: set region failed. Check that REGION describes a valid range" << endl;
 111                     reader.Close();
 112                     return 1;
 113                 }
 114
 115                 // everything checks out, just iterate through specified region, counting alignments
 116                 while ( reader.GetNextAlignmentCore(al) )
 117                     ++alignmentCount;
 118             }
 119
 120             // no index data available, we have to iterate through until we
 121             // find overlapping alignments
 122             else {
 123                 while ( reader.GetNextAlignmentCore(al) ) {
 124                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 125                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 126                     {
 127                         ++alignmentCount;
 128                     }
 129                 }
 130             }
 131         }
 132
 133         // error parsing REGION string
 134         else {
 135             cerr << "bamtools count ERROR: could not parse REGION - " << m_settings->Region << endl;
 136             cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 137                  << endl;
 138             reader.Close();
 139             return 1;
 140         }
 141     }
 142
 143     // print results
 144     cout << alignmentCount << endl;
 145
 146     // clean up & exit
 147     reader.Close();
 148     return 0;
 149 }