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