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     if (!reader.Open(m_settings->InputFiles, false, true)) {
  84         cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
  85         return 1;
  86     }
  87
  88     // alignment counter
  89     BamAlignment al;
  90     int alignmentCount(0);
  91
  92     // if no region specified, count entire file
  93     if ( !m_settings->HasRegion ) {
  94         while ( reader.GetNextAlignmentCore(al) )
  95             ++alignmentCount;
  96     }
  97
  98     // otherwise attempt to use region as constraint
  99     else {
 100
 101         // if region string parses OK
 102         BamRegion region;
 103         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 104
 105             // attempt to re-open reader with index files
 106             reader.Close();
 107             bool openedOK = reader.Open(m_settings->InputFiles, true, true );
 108
 109             // if error
 110             if ( !openedOK ) {
 111                 cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
 112                 return 1;
 113             }
 114
 115             // if index data available, we can use SetRegion
 116             if ( reader.IsIndexLoaded() ) {
 117
 118                 // attempt to use SetRegion(), if failed report error
 119                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 120                     cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
 121                     reader.Close();
 122                     return 1;
 123                 }
 124
 125                 // everything checks out, just iterate through specified region, counting alignments
 126                 while ( reader.GetNextAlignmentCore(al) )
 127                     ++alignmentCount;
 128             }
 129
 130             // no index data available, we have to iterate through until we
 131             // find overlapping alignments
 132             else {
 133                 while( reader.GetNextAlignmentCore(al) ) {
 134                     if ( (al.RefID >= region.LeftRefID)  && ( (al.Position + al.Length) >= region.LeftPosition ) &&
 135                           (al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
 136                     {
 137                         ++alignmentCount;
 138                     }
 139                 }
 140             }
 141         }
 142
 143         // error parsing REGION string
 144         else {
 145             cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
 146             cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 147             reader.Close();
 148             return 1;
 149         }
 150     }
 151
 152     // print results
 153     cout << alignmentCount << endl;
 154
 155     // clean & exit
 156     reader.Close();
 157     return 0;
 158 }