bamtools_utilities.cpp

   1 // ***************************************************************************
   2 // bamtools_utilities.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 // Provides general utilities used by BamTools sub-tools.
   9 // ***************************************************************************
  10
  11 #include <cstdlib>
  12 #include "bamtools_utilities.h"
  13 #include "BamReader.h"
  14 #include "BamMultiReader.h"
  15
  16 using namespace std;
  17 using namespace BamTools;
  18
  19 // Parses a region string, does validation (valid ID's, positions), stores in Region struct
  20 // Returns success (true/false)
  21 bool Utilities::ParseRegionString(const std::string& regionString, const BamReader& reader, Region& region) {
  22
  23     // -------------------------------
  24     // parse region string
  25
  26     // check first for empty string
  27     if ( regionString.empty() )
  28         return false;
  29
  30     // non-empty string, look for a colom
  31     size_t foundFirstColon = regionString.find(':');
  32
  33     // store chrom strings, and numeric positions
  34     string startChrom;
  35     string stopChrom;
  36     int startPos;
  37     int stopPos;
  38
  39     // no colon found
  40     // going to use entire contents of requested chromosome
  41     // just store entire region string as startChrom name
  42     // use BamReader methods to check if its valid for current BAM file
  43     if ( foundFirstColon == string::npos ) {
  44         startChrom = regionString;
  45         startPos   = 0;
  46         stopChrom  = regionString;
  47         stopPos    = -1;
  48     }
  49
  50     // colon found, so we at least have some sort of startPos requested
  51     else {
  52
  53         // store start chrom from beginning to first colon
  54         startChrom = regionString.substr(0,foundFirstColon);
  55
  56         // look for ".." after the colon
  57         size_t foundRangeDots = regionString.find("..", foundFirstColon+1);
  58
  59         // no dots found
  60         // so we have a startPos but no range
  61         // store contents before colon as startChrom, after as startPos
  62         if ( foundRangeDots == string::npos ) {
  63             startPos   = atoi( regionString.substr(foundFirstColon+1).c_str() );
  64             stopChrom  = startChrom;
  65             stopPos    = -1;
  66         }
  67
  68         // ".." found, so we have some sort of range selected
  69         else {
  70
  71             // store startPos between first colon and range dots ".."
  72             startPos = atoi( regionString.substr(foundFirstColon+1, foundRangeDots-foundFirstColon-1).c_str() );
  73
  74             // look for second colon
  75             size_t foundSecondColon = regionString.find(':', foundRangeDots+1);
  76
  77             // no second colon found
  78             // so we have a "standard" chrom:start..stop input format (on single chrom)
  79             if ( foundSecondColon == string::npos ) {
  80                 stopChrom  = startChrom;
  81                 stopPos    = atoi( regionString.substr(foundRangeDots+2).c_str() );
  82             }
  83
  84             // second colon found
  85             // so we have a range requested across 2 chrom's
  86             else {
  87                 stopChrom  = regionString.substr(foundRangeDots+2, foundSecondColon-(foundRangeDots+2));
  88                 stopPos    = atoi( regionString.substr(foundSecondColon+1).c_str() );
  89             }
  90         }
  91     }
  92
  93     // -------------------------------
  94     // validate reference IDs & genomic positions
  95
  96     const RefVector references = reader.GetReferenceData();
  97
  98     // if startRefID not found, return false
  99     int startRefID = reader.GetReferenceID(startChrom);
 100     if ( startRefID == (int)references.size() ) return false;
 101
 102     // if startPos is larger than reference, return false
 103     const RefData& startReference = references.at(startRefID);
 104     if ( startPos > startReference.RefLength ) return false;
 105
 106     // if stopRefID not found, return false
 107     int stopRefID = reader.GetReferenceID(stopChrom);
 108     if ( stopRefID == (int)references.size() ) return false;
 109
 110     // if stopPosition larger than reference, return false
 111     const RefData& stopReference = references.at(stopRefID);
 112     if ( stopPos > stopReference.RefLength ) return false;
 113
 114     // if no stopPosition specified, set to reference end
 115     if ( stopPos == -1 ) stopPos = stopReference.RefLength;
 116
 117     // -------------------------------
 118     // set up Region struct & return
 119
 120     region.StartChromID  = startRefID;
 121     region.StopChromID   = stopRefID;
 122     region.StartPosition = startPos;
 123     region.StopPosition  = stopPos;
 124     return true;
 125 }
 126
 127 // Same as ParseRegionString() above, but accepts a BamMultiReader
 128 bool Utilities::ParseRegionString(const std::string& regionString, const BamMultiReader& reader, Region& region) {
 129
 130     // -------------------------------
 131     // parse region string
 132
 133     // check first for empty string
 134     if ( regionString.empty() )
 135         return false;
 136
 137     // non-empty string, look for a colom
 138     size_t foundFirstColon = regionString.find(':');
 139
 140     // store chrom strings, and numeric positions
 141     string startChrom;
 142     string stopChrom;
 143     int startPos;
 144     int stopPos;
 145
 146     // no colon found
 147     // going to use entire contents of requested chromosome
 148     // just store entire region string as startChrom name
 149     // use BamReader methods to check if its valid for current BAM file
 150     if ( foundFirstColon == string::npos ) {
 151         startChrom = regionString;
 152         startPos   = 0;
 153         stopChrom  = regionString;
 154         stopPos    = -1;
 155     }
 156
 157     // colon found, so we at least have some sort of startPos requested
 158     else {
 159
 160         // store start chrom from beginning to first colon
 161         startChrom = regionString.substr(0,foundFirstColon);
 162
 163         // look for ".." after the colon
 164         size_t foundRangeDots = regionString.find("..", foundFirstColon+1);
 165
 166         // no dots found
 167         // so we have a startPos but no range
 168         // store contents before colon as startChrom, after as startPos
 169         if ( foundRangeDots == string::npos ) {
 170             startPos   = atoi( regionString.substr(foundFirstColon+1).c_str() );
 171             stopChrom  = startChrom;
 172             stopPos    = -1;
 173         }
 174
 175         // ".." found, so we have some sort of range selected
 176         else {
 177
 178             // store startPos between first colon and range dots ".."
 179             startPos = atoi( regionString.substr(foundFirstColon+1, foundRangeDots-foundFirstColon-1).c_str() );
 180
 181             // look for second colon
 182             size_t foundSecondColon = regionString.find(':', foundRangeDots+1);
 183
 184             // no second colon found
 185             // so we have a "standard" chrom:start..stop input format (on single chrom)
 186             if ( foundSecondColon == string::npos ) {
 187                 stopChrom  = startChrom;
 188                 stopPos    = atoi( regionString.substr(foundRangeDots+2).c_str() );
 189             }
 190
 191             // second colon found
 192             // so we have a range requested across 2 chrom's
 193             else {
 194                 stopChrom  = regionString.substr(foundRangeDots+2, foundSecondColon-(foundRangeDots+2));
 195                 stopPos    = atoi( regionString.substr(foundSecondColon+1).c_str() );
 196             }
 197         }
 198     }
 199
 200     // -------------------------------
 201     // validate reference IDs & genomic positions
 202
 203     const RefVector references = reader.GetReferenceData();
 204
 205     // if startRefID not found, return false
 206     int startRefID = reader.GetReferenceID(startChrom);
 207     if ( startRefID == (int)references.size() ) return false;
 208
 209     // if startPos is larger than reference, return false
 210     const RefData& startReference = references.at(startRefID);
 211     if ( startPos > startReference.RefLength ) return false;
 212
 213     // if stopRefID not found, return false
 214     int stopRefID = reader.GetReferenceID(stopChrom);
 215     if ( stopRefID == (int)references.size() ) return false;
 216
 217     // if stopPosition larger than reference, return false
 218     const RefData& stopReference = references.at(stopRefID);
 219     if ( stopPos > stopReference.RefLength ) return false;
 220
 221     // if no stopPosition specified, set to reference end
 222     if ( stopPos == -1 ) stopPos = stopReference.RefLength;
 223
 224     // -------------------------------
 225     // set up Region struct & return
 226
 227     region.StartChromID  = startRefID;
 228     region.StopChromID   = stopRefID;
 229     region.StartPosition = startPos;
 230     region.StopPosition  = stopPos;
 231     return true;
 232 }