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: 2 June 2010
   7 // ---------------------------------------------------------------------------
   8 // Prints alignment count for BAM file
   9 //
  10 // ** Expand to multiple??
  11 //
  12 // ***************************************************************************
  13
  14 #include <iostream>
  15 #include <sstream>
  16 #include <string>
  17 #include <vector>
  18
  19 #include "bamtools_count.h"
  20 #include "bamtools_options.h"
  21 #include "bamtools_utilities.h"
  22 #include "BamReader.h"
  23
  24 using namespace std;
  25 using namespace BamTools;
  26
  27 // ---------------------------------------------
  28 // CountSettings implementation
  29
  30 struct CountTool::CountSettings {
  31
  32     // flags
  33     bool HasBamIndexFilename;
  34     bool HasInputBamFilename;
  35     bool HasRegion;
  36
  37     // filenames
  38     string BamIndexFilename;
  39     string InputBamFilename;
  40     string Region;
  41
  42     // constructor
  43     CountSettings(void)
  44         : HasBamIndexFilename(false)
  45         , HasInputBamFilename(false)
  46         , HasRegion(false)
  47         , InputBamFilename(Options::StandardIn())
  48     { }
  49 };
  50
  51 // ---------------------------------------------
  52 // CountTool implementation
  53
  54 CountTool::CountTool(void)
  55     : AbstractTool()
  56     , m_settings(new CountSettings)
  57 {
  58     // set program details
  59     Options::SetProgramInfo("bamtools count", "prints alignment counts for a BAM file", "-in <filename> [-region <REGION> [-index <filename>]]");
  60
  61     // set up options
  62     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  63     Options::AddValueOption("-in",    "BAM filename",        "the input BAM file", "", m_settings->HasInputBamFilename, m_settings->InputBamFilename, IO_Opts, Options::StandardIn());
  64     Options::AddValueOption("-index", "BAM index filename", "the BAM index file",  "", m_settings->HasBamIndexFilename, m_settings->BamIndexFilename, IO_Opts);
  65
  66     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  67     Options::AddValueOption("-region", "REGION", "genomic region. Index file is recommended for better performance. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
  68 }
  69
  70 CountTool::~CountTool(void) {
  71     delete m_settings;
  72     m_settings = 0;
  73 }
  74
  75 int CountTool::Help(void) {
  76     Options::DisplayHelp();
  77     return 0;
  78 }
  79
  80 int CountTool::Run(int argc, char* argv[]) {
  81
  82     // parse command line arguments
  83     Options::Parse(argc, argv, 1);
  84
  85     //open our BAM reader
  86     BamReader reader;
  87     reader.Open(m_settings->InputBamFilename);
  88
  89     // alignment counter
  90     int alignmentCount(0);
  91
  92     // set up error handling
  93     ostringstream errorStream("");
  94     bool foundError(false);
  95
  96     // if no region specified, count entire file
  97     if ( !m_settings->HasRegion ) {
  98         BamAlignment al;
  99         while ( reader.GetNextAlignment(al) )
 100             ++alignmentCount;
 101     }
 102
 103     // more complicated - region specified
 104     else {
 105
 106         Region region;
 107         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 108
 109             // has index, we can jump directly to
 110             if ( m_settings->HasBamIndexFilename ) {
 111
 112                 // this is kind of a hack...?
 113                 // re-open BamReader, this time with the index file
 114                 reader.Close();
 115                 reader.Open(m_settings->InputBamFilename, m_settings->BamIndexFilename);
 116
 117                 // attempt Jump(), catch error
 118                 if ( !reader.Jump(region.StartChromID, region.StartPosition) ) {
 119                     foundError = true;
 120                     errorStream << "Could not jump to desired REGION: " << m_settings->Region << endl;
 121                 }
 122             }
 123
 124             else {
 125
 126                 // read through sequentially, until first overlapping read is found
 127                 BamAlignment al;
 128                 bool alignmentFound(false);
 129                 while( reader.GetNextAlignment(al) ) {
 130                     if ( (al.RefID == region.StartChromID ) && ( (al.Position + al.Length) >= region.StartPosition) ) {
 131                         alignmentFound = true;
 132                         break;
 133                     }
 134                 }
 135
 136                 // if overlapping alignment found (not EOF), increment counter
 137                 if ( alignmentFound) ++ alignmentCount;
 138             }
 139
 140             // -----------------------------
 141             // count alignments until stop hit
 142
 143             if ( !foundError ) {
 144                 // while valid alignment AND
 145                 // ( either current ref is before stopRefID OR
 146                 //   current ref stopRefID but we're before stopPos )
 147                 BamAlignment al;
 148                 while ( reader.GetNextAlignment(al) && ((al.RefID < region.StopChromID ) || (al.RefID == region.StopChromID && al.Position <= region.StopPosition)) )
 149                     ++alignmentCount;
 150             }
 151
 152         } else {
 153             foundError = true;
 154             errorStream << "Could not parse REGION: " << m_settings->Region << endl;
 155             errorStream << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 156         }
 157     }
 158
 159     // print errors OR results
 160     if ( foundError )
 161         cerr << errorStream.str() << endl;
 162     else
 163         cout << alignmentCount << endl;
 164
 165     // clean & exit
 166     reader.Close();
 167     return (int)foundError;
 168 }