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 alignment counts for a BAM file", "[-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
  58     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  59     Options::AddValueOption("-region", "REGION", "genomic region. Index file is required and is read automatically if it exists as <filename>.bai or <filename>.bti. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
  60 }
  61
  62 CountTool::~CountTool(void) {
  63     delete m_settings;
  64     m_settings = 0;
  65 }
  66
  67 int CountTool::Help(void) {
  68     Options::DisplayHelp();
  69     return 0;
  70 }
  71
  72 int CountTool::Run(int argc, char* argv[]) {
  73
  74     // parse command line arguments
  75     Options::Parse(argc, argv, 1);
  76
  77     // if no '-in' args supplied, default to stdin
  78     if ( !m_settings->HasInput )
  79         m_settings->InputFiles.push_back(Options::StandardIn());
  80
  81     // open reader without index
  82     BamMultiReader reader;
  83     reader.Open(m_settings->InputFiles, false, true);
  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 re-open reader with index files
 103             reader.Close();
 104             bool openedOK = reader.Open(m_settings->InputFiles, true, true );
 105
 106             // if error
 107             if ( !openedOK ) {
 108                 cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
 109                 return 1;
 110             }
 111
 112             // if index data available, we can use SetRegion
 113             if ( reader.IsIndexLoaded() ) {
 114
 115                 // attempt to use SetRegion(), if failed report error
 116                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 117                     cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
 118                     reader.Close();
 119                     return 1;
 120                 }
 121
 122                 // everything checks out, just iterate through specified region, counting alignments
 123                 while ( reader.GetNextAlignmentCore(al) )
 124                     ++alignmentCount;
 125             }
 126
 127             // no index data available, we have to iterate through until we
 128             // find overlapping alignments
 129             else {
 130                 while( reader.GetNextAlignmentCore(al) ) {
 131                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 132                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 133                     {
 134                         ++alignmentCount;
 135                     }
 136                 }
 137             }
 138         }
 139
 140         // error parsing REGION string
 141         else {
 142             cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
 143             cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 144             reader.Close();
 145             return 1;
 146         }
 147     }
 148
 149     // print results
 150     cout << alignmentCount << endl;
 151
 152     // clean & exit
 153     reader.Close();
 154     return 0;
 155 }