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         BamRegion 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                 if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
 124                    foundError = true;
 125                    errorStream << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
 126                 }
 127
 128             }
 129
 130             else {
 131
 132                 // read through sequentially, until first overlapping read is found
 133                 BamAlignment al;
 134                 bool alignmentFound(false);
 135                 while( reader.GetNextAlignment(al) ) {
 136                     if ( (al.RefID == region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) ) {
 137                         alignmentFound = true;
 138                         break;
 139                     }
 140                 }
 141
 142                 // if overlapping alignment found (not EOF), increment counter
 143                 if ( alignmentFound ) ++ alignmentCount;
 144             }
 145
 146             // -----------------------------
 147             // count alignments until stop hit
 148
 149             if ( !foundError ) {
 150                 // while valid alignment AND
 151                 // ( either current ref is before stopRefID OR
 152                 //   current ref stopRefID but we're before stopPos )
 153                 BamAlignment al;
 154                 //while ( reader.GetNextAlignment(al) && ((al.RefID < region.StopChromID ) || (al.RefID == region.StopChromID && al.Position <= region.StopPosition)) )
 155                 //    ++alignmentCount;
 156                 //
 157                 while ( reader.GetNextAlignment(al) )
 158                     ++alignmentCount;
 159
 160             }
 161
 162         } else {
 163             foundError = true;
 164             errorStream << "Could not parse REGION: " << m_settings->Region << endl;
 165             errorStream << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 166         }
 167     }
 168
 169     // print errors OR results
 170     if ( foundError )
 171         cerr << errorStream.str() << endl;
 172     else
 173         cout << alignmentCount << endl;
 174
 175     // clean & exit
 176     reader.Close();
 177     return (int)foundError;
 178 }