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