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 optimal performance.", "", 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     } else {
 102
 103         // parse region string for desired region
 104         string startChrom;
 105         string stopChrom;
 106         int startPos;
 107         int stopPos;
 108         if ( Utilities::ParseRegionString(m_settings->Region, startChrom, startPos, stopChrom, stopPos) ) {
 109
 110             const RefVector references = reader.GetReferenceData();
 111
 112             // -------------------------------
 113             // validate start ref & position
 114
 115             int startRefID = reader.GetReferenceID(startChrom);
 116             cout << "startRefID: " << startRefID << endl;
 117
 118             // startRefID not found
 119             if ( startRefID == (int)references.size() ) {
 120                 foundError = true;
 121                 errorStream << "Start chrom: " << startChrom << " not found in BAM file." << endl;
 122             }
 123
 124             // valid startRefID, check startPos
 125             else {
 126                 const RefData& reference = references.at(startRefID);
 127
 128                 // startPos too large
 129                 if ( startPos > reference.RefLength ) {
 130                     foundError = true;
 131                     errorStream << "Start pos: " << startPos << " is larger than expected." << endl;
 132                 }
 133             }
 134
 135             // -------------------------------
 136             // validate stop ref & position
 137
 138             int stopRefID = reader.GetReferenceID(stopChrom);
 139
 140             // skip if error already found
 141             if ( !foundError ) {
 142
 143                 // stopRefID not found
 144                 if ( stopRefID == (int)references.size() ) {
 145                     foundError = true;
 146                     errorStream << "Stop chrom: " << stopChrom << " not found in BAM file." << endl;
 147                 }
 148
 149                 // valid stopRefID, check stopPos
 150                 else {
 151                     const RefData& reference = references.at(stopRefID);
 152
 153                     // stopPos too large
 154                     if ( stopPos > reference.RefLength ) {
 155                         foundError = true;
 156                         errorStream << "Stop pos: " << stopPos << " is larger than expected." << endl;
 157                     }
 158
 159                     // no stopPos specified, set to reference end
 160                     else if ( stopPos == -1 ) {
 161                         stopPos = reference.RefLength;
 162                     }
 163                 }
 164             }
 165
 166             // -------------------------------
 167             // if refs & positions valid, retrieve data
 168
 169             if ( !foundError ) {
 170
 171                 // has index, we can jump directly to
 172                 if ( m_settings->HasBamIndexFilename ) {
 173
 174                     // this is kind of a hack...?
 175                     // re-open BamReader, this time with the index file
 176                     reader.Close();
 177                     reader.Open(m_settings->InputBamFilename, m_settings->BamIndexFilename);
 178
 179                     // attempt Jump(), catch error
 180                     if ( !reader.Jump(startRefID, startPos) ) {
 181                         foundError = true;
 182                         errorStream << "Could not jump to desired REGION: " << m_settings->Region << endl;
 183                     }
 184                 }
 185
 186                 else {
 187
 188                     // read through sequentially, until first overlapping read is found
 189                     BamAlignment al;
 190                     bool alignmentFound(false);
 191                     while( reader.GetNextAlignment(al) ) {
 192                         if ( (al.RefID == startRefID ) && ( (al.Position + al.Length) >= startPos) ) {
 193                             alignmentFound = true;
 194                             break;
 195                         }
 196                     }
 197
 198                     // if overlapping alignment found (not EOF), increment counter
 199                     if ( alignmentFound) ++ alignmentCount;
 200                 }
 201
 202                 // -----------------------------
 203                 // count alignments until
 204
 205                 if ( !foundError ) {
 206                     // while valid alignment AND
 207                     // ( either current ref is before stopRefID OR
 208                     //   current ref stopRefID but we're before stopPos )
 209                     BamAlignment al;
 210                     while ( reader.GetNextAlignment(al) && ((al.RefID < stopRefID ) || (al.RefID == stopRefID && al.Position <= stopPos)) )
 211                         ++alignmentCount;
 212                 }
 213             }
 214         }
 215
 216         // could not parse REGION string, set error
 217         else {
 218             foundError = true;
 219             errorStream << "Could not parse desired REGION: " << m_settings->Region << endl;
 220             errorStream << "Please see README for details on accepted REGION formats" << endl;
 221         }
 222     }
 223
 224     // print errors OR results
 225     if ( foundError )
 226         cerr << errorStream.str() << endl;
 227     else
 228         cout << alignmentCount << endl;
 229
 230     // clean & exit
 231     reader.Close();
 232     return (int)foundError;
 233 }