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: 3 September 2010
   7 // ---------------------------------------------------------------------------
   8 // Prints alignment count for BAM file(s)
   9 // ***************************************************************************
  10
  11 #include <iostream>
  12 #include <string>
  13 #include <vector>
  14
  15 #include "bamtools_count.h"
  16 #include "bamtools_options.h"
  17 #include "bamtools_utilities.h"
  18 #include "BamReader.h"
  19 #include "BamMultiReader.h"
  20
  21 using namespace std;
  22 using namespace BamTools;
  23
  24 // ---------------------------------------------
  25 // CountSettings implementation
  26
  27 struct CountTool::CountSettings {
  28
  29     // flags
  30     bool HasInput;
  31     bool HasRegion;
  32
  33     // filenames
  34     vector<string> InputFiles;
  35     string Region;
  36
  37     // constructor
  38     CountSettings(void)
  39         : HasInput(false)
  40         , HasRegion(false)
  41     { }
  42 };
  43
  44 // ---------------------------------------------
  45 // CountTool implementation
  46
  47 CountTool::CountTool(void)
  48     : AbstractTool()
  49     , m_settings(new CountSettings)
  50 {
  51     // set program details
  52     Options::SetProgramInfo("bamtools count", "prints number of alignments in BAM file(s)", "[-in <filename> -in <filename> ...] [-region <REGION>]");
  53
  54     // set up options
  55     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  56     Options::AddValueOption("-in",     "BAM filename", "the input BAM file(s)", "", m_settings->HasInput,  m_settings->InputFiles, IO_Opts, Options::StandardIn());
  57     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);
  58 }
  59
  60 CountTool::~CountTool(void) {
  61     delete m_settings;
  62     m_settings = 0;
  63 }
  64
  65 int CountTool::Help(void) {
  66     Options::DisplayHelp();
  67     return 0;
  68 }
  69
  70 int CountTool::Run(int argc, char* argv[]) {
  71
  72     // parse command line arguments
  73     Options::Parse(argc, argv, 1);
  74
  75     // if no '-in' args supplied, default to stdin
  76     if ( !m_settings->HasInput )
  77         m_settings->InputFiles.push_back(Options::StandardIn());
  78
  79     // open reader without index
  80     BamMultiReader reader;
  81     if (!reader.Open(m_settings->InputFiles, false, true)) {
  82         cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
  83         return 1;
  84     }
  85
  86     // alignment counter
  87     BamAlignment al;
  88     int alignmentCount(0);
  89
  90     // if no region specified, count entire file
  91     if ( !m_settings->HasRegion ) {
  92         while ( reader.GetNextAlignmentCore(al) )
  93             ++alignmentCount;
  94     }
  95
  96     // otherwise attempt to use region as constraint
  97     else {
  98
  99         // if region string parses OK
 100         BamRegion region;
 101         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 102
 103             // attempt to re-open reader with index files
 104             reader.Close();
 105             bool openedOK = reader.Open(m_settings->InputFiles, true, true );
 106
 107             // if error
 108             if ( !openedOK ) {
 109                 cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
 110                 return 1;
 111             }
 112
 113             // if index data available, we can use SetRegion
 114             if ( reader.IsIndexLoaded() ) {
 115
 116                 // attempt to use SetRegion(), if failed report error
 117                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 118                     cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
 119                     reader.Close();
 120                     return 1;
 121                 }
 122
 123                 // everything checks out, just iterate through specified region, counting alignments
 124                 while ( reader.GetNextAlignmentCore(al) )
 125                     ++alignmentCount;
 126             }
 127
 128             // no index data available, we have to iterate through until we
 129             // find overlapping alignments
 130             else {
 131                 while( reader.GetNextAlignmentCore(al) ) {
 132                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 133                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 134                     {
 135                         ++alignmentCount;
 136                     }
 137                 }
 138             }
 139         }
 140
 141         // error parsing REGION string
 142         else {
 143             cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
 144             cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 145             reader.Close();
 146             return 1;
 147         }
 148     }
 149
 150     // print results
 151     cout << alignmentCount << endl;
 152
 153     // clean & exit
 154     reader.Close();
 155     return 0;
 156 }