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 // ---------------------------------------------------------------------------
   5 // Last modified: 9 June 2011
   6 // ---------------------------------------------------------------------------
   7 // Provides general utilities used by BamTools sub-tools.
   8 // ***************************************************************************
   9
  10 #include <api/BamMultiReader.h>
  11 #include <api/BamReader.h>
  12 #include <utils/bamtools_utilities.h>
  13 using namespace BamTools;
  14
  15 #include <algorithm>
  16 #include <cstdlib>
  17 #include <cstring>
  18 #include <fstream>
  19 #include <iostream>
  20 #include <sstream>
  21 using namespace std;
  22
  23 namespace BamTools {
  24
  25 const char REVCOMP_LOOKUP[] = {'T',  0,  'G', 'H',
  26                                 0,   0,  'C', 'D',
  27                                 0,   0,   0,   0,
  28                                'K', 'N',  0,   0,
  29                                 0,  'Y', 'W', 'A',
  30                                'A', 'B', 'S', 'X',
  31                                'R',  0 };
  32
  33 } // namespace BamTools
  34
  35 // returns true if 'source' contains 'pattern'
  36 bool Utilities::Contains(const string& source, const string& pattern) {
  37     return ( source.find(pattern) != string::npos );
  38 }
  39
  40 // returns true if 'source' contains 'c'
  41 bool Utilities::Contains(const std::string &source, const char c) {
  42     return ( source.find(c) != string::npos );
  43 }
  44
  45 // returns true if 'source' ends with 'pattern'
  46 bool Utilities::EndsWith(const string& source, const string& pattern) {
  47     return ( source.find(pattern) == (source.length() - pattern.length()) );
  48 }
  49
  50 // returns true if 'source' ends with 'c'
  51 bool Utilities::EndsWith(const std::string& source, const char c) {
  52     return ( source.find(c) == (source.length() - 1) );
  53 }
  54
  55 // check if a file exists
  56 bool Utilities::FileExists(const string& filename) {
  57     ifstream f(filename.c_str(), ifstream::in);
  58     return !f.fail();
  59 }
  60
  61 // Parses a region string, does validation (valid ID's, positions), stores in Region struct
  62 // Returns success (true/false)
  63 bool Utilities::ParseRegionString(const string& regionString,
  64                                   const BamReader& reader,
  65                                   BamRegion& region)
  66 {
  67     // -------------------------------
  68     // parse region string
  69
  70     // check first for empty string
  71     if ( regionString.empty() )
  72         return false;
  73
  74     // non-empty string, look for a colom
  75     size_t foundFirstColon = regionString.find(':');
  76
  77     // store chrom strings, and numeric positions
  78     string startChrom;
  79     string stopChrom;
  80     int startPos;
  81     int stopPos;
  82
  83     // no colon found
  84     // going to use entire contents of requested chromosome
  85     // just store entire region string as startChrom name
  86     // use BamReader methods to check if its valid for current BAM file
  87     if ( foundFirstColon == string::npos ) {
  88         startChrom = regionString;
  89         startPos   = 0;
  90         stopChrom  = regionString;
  91         stopPos    = -1;
  92     }
  93
  94     // colon found, so we at least have some sort of startPos requested
  95     else {
  96
  97         // store start chrom from beginning to first colon
  98         startChrom = regionString.substr(0,foundFirstColon);
  99
 100         // look for ".." after the colon
 101         size_t foundRangeDots = regionString.find("..", foundFirstColon+1);
 102
 103         // no dots found
 104         // so we have a startPos but no range
 105         // store contents before colon as startChrom, after as startPos
 106         if ( foundRangeDots == string::npos ) {
 107             startPos   = atoi( regionString.substr(foundFirstColon+1).c_str() );
 108             stopChrom  = startChrom;
 109             stopPos    = -1;
 110         }
 111
 112         // ".." found, so we have some sort of range selected
 113         else {
 114
 115             // store startPos between first colon and range dots ".."
 116             startPos = atoi( regionString.substr(foundFirstColon+1, foundRangeDots-foundFirstColon-1).c_str() );
 117
 118             // look for second colon
 119             size_t foundSecondColon = regionString.find(':', foundRangeDots+1);
 120
 121             // no second colon found
 122             // so we have a "standard" chrom:start..stop input format (on single chrom)
 123             if ( foundSecondColon == string::npos ) {
 124                 stopChrom  = startChrom;
 125                 stopPos    = atoi( regionString.substr(foundRangeDots+2).c_str() );
 126             }
 127
 128             // second colon found
 129             // so we have a range requested across 2 chrom's
 130             else {
 131                 stopChrom  = regionString.substr(foundRangeDots+2, foundSecondColon-(foundRangeDots+2));
 132                 stopPos    = atoi( regionString.substr(foundSecondColon+1).c_str() );
 133             }
 134         }
 135     }
 136
 137     // -------------------------------
 138     // validate reference IDs & genomic positions
 139
 140     const RefVector references = reader.GetReferenceData();
 141
 142     // if startRefID not found, return false
 143     int startRefID = reader.GetReferenceID(startChrom);
 144     if ( startRefID == (int)references.size() ) return false;
 145
 146     // if startPos is larger than reference, return false
 147     const RefData& startReference = references.at(startRefID);
 148     if ( startPos > startReference.RefLength ) return false;
 149
 150     // if stopRefID not found, return false
 151     int stopRefID = reader.GetReferenceID(stopChrom);
 152     if ( stopRefID == (int)references.size() ) return false;
 153
 154     // if stopPosition larger than reference, return false
 155     const RefData& stopReference = references.at(stopRefID);
 156     if ( stopPos > stopReference.RefLength ) return false;
 157
 158     // if no stopPosition specified, set to reference end
 159     if ( stopPos == -1 ) stopPos = stopReference.RefLength;
 160
 161     // -------------------------------
 162     // set up Region struct & return
 163
 164     region.LeftRefID = startRefID;
 165     region.LeftPosition = startPos;
 166     region.RightRefID = stopRefID;;
 167     region.RightPosition = stopPos;
 168     return true;
 169 }
 170
 171 // Same as ParseRegionString() above, but accepts a BamMultiReader
 172 bool Utilities::ParseRegionString(const string& regionString,
 173                                   const BamMultiReader& reader,
 174                                   BamRegion& region)
 175 {
 176     // -------------------------------
 177     // parse region string
 178
 179     // check first for empty string
 180     if ( regionString.empty() )
 181         return false;
 182
 183     // non-empty string, look for a colom
 184     size_t foundFirstColon = regionString.find(':');
 185
 186     // store chrom strings, and numeric positions
 187     string startChrom;
 188     string stopChrom;
 189     int startPos;
 190     int stopPos;
 191
 192     // no colon found
 193     // going to use entire contents of requested chromosome
 194     // just store entire region string as startChrom name
 195     // use BamReader methods to check if its valid for current BAM file
 196     if ( foundFirstColon == string::npos ) {
 197         startChrom = regionString;
 198         startPos   = 0;
 199         stopChrom  = regionString;
 200         stopPos    = -1;
 201     }
 202
 203     // colon found, so we at least have some sort of startPos requested
 204     else {
 205
 206         // store start chrom from beginning to first colon
 207         startChrom = regionString.substr(0,foundFirstColon);
 208
 209         // look for ".." after the colon
 210         size_t foundRangeDots = regionString.find("..", foundFirstColon+1);
 211
 212         // no dots found
 213         // so we have a startPos but no range
 214         // store contents before colon as startChrom, after as startPos
 215         if ( foundRangeDots == string::npos ) {
 216             startPos   = atoi( regionString.substr(foundFirstColon+1).c_str() );
 217             stopChrom  = startChrom;
 218             stopPos    = -1;
 219         }
 220
 221         // ".." found, so we have some sort of range selected
 222         else {
 223
 224             // store startPos between first colon and range dots ".."
 225             startPos = atoi( regionString.substr(foundFirstColon+1, foundRangeDots-foundFirstColon-1).c_str() );
 226
 227             // look for second colon
 228             size_t foundSecondColon = regionString.find(':', foundRangeDots+1);
 229
 230             // no second colon found
 231             // so we have a "standard" chrom:start..stop input format (on single chrom)
 232             if ( foundSecondColon == string::npos ) {
 233                 stopChrom  = startChrom;
 234                 stopPos    = atoi( regionString.substr(foundRangeDots+2).c_str() );
 235             }
 236
 237             // second colon found
 238             // so we have a range requested across 2 chrom's
 239             else {
 240                 stopChrom  = regionString.substr(foundRangeDots+2, foundSecondColon-(foundRangeDots+2));
 241                 stopPos    = atoi( regionString.substr(foundSecondColon+1).c_str() );
 242             }
 243         }
 244     }
 245
 246     // -------------------------------
 247     // validate reference IDs & genomic positions
 248
 249     const RefVector references = reader.GetReferenceData();
 250
 251     // if startRefID not found, return false
 252     int startRefID = reader.GetReferenceID(startChrom);
 253     if ( startRefID == (int)references.size() ) return false;
 254
 255     // if startPos is larger than reference, return false
 256     const RefData& startReference = references.at(startRefID);
 257     if ( startPos > startReference.RefLength ) return false;
 258
 259     // if stopRefID not found, return false
 260     int stopRefID = reader.GetReferenceID(stopChrom);
 261     if ( stopRefID == (int)references.size() ) return false;
 262
 263     // if stopPosition larger than reference, return false
 264     const RefData& stopReference = references.at(stopRefID);
 265     if ( stopPos > stopReference.RefLength ) return false;
 266
 267     // if no stopPosition specified, set to reference end
 268     if ( stopPos == -1 ) stopPos = stopReference.RefLength;
 269
 270     // -------------------------------
 271     // set up Region struct & return
 272
 273     region.LeftRefID = startRefID;
 274     region.LeftPosition = startPos;
 275     region.RightRefID = stopRefID;;
 276     region.RightPosition = stopPos;
 277     return true;
 278 }
 279
 280 void Utilities::Reverse(string& sequence) {
 281     reverse(sequence.begin(), sequence.end());
 282 }
 283
 284 void Utilities::ReverseComplement(string& sequence) {
 285
 286     // do complement, in-place
 287     size_t seqLength = sequence.length();
 288     for ( size_t i = 0; i < seqLength; ++i )
 289         sequence.replace(i, 1, 1, REVCOMP_LOOKUP[(int)sequence.at(i) - 65]);
 290
 291     // reverse it
 292     Reverse(sequence);
 293 }
 294
 295 vector<string> Utilities::Split(const string& source, const char delim) {
 296
 297     stringstream ss(source);
 298     string field;
 299     vector<string> fields;
 300
 301     while ( getline(ss, field, delim) )
 302         fields.push_back(field);
 303     return fields;
 304 }
 305
 306 vector<string> Utilities::Split(const string& source, const string& delims) {
 307
 308     vector<string> fields;
 309
 310     char* tok;
 311     char* cchars = new char[source.size()+1];
 312     char* cstr = &cchars[0];
 313     strcpy(cstr, source.c_str());
 314     tok = strtok(cstr, delims.c_str());
 315     while (tok != NULL) {
 316         fields.push_back(tok);
 317         tok = strtok(NULL, delims.c_str());
 318     }
 319
 320     delete[] cchars;
 321
 322     return fields;
 323 }
 324
 325 // returns true if 'source' starts with 'pattern'
 326 bool Utilities::StartsWith(const string& source, const string& pattern) {
 327     return ( source.find(pattern) == 0 );
 328 }
 329
 330 // returns true if 'source' starts with 'c'
 331 bool Utilities::StartsWith(const std::string &source, const char c) {
 332     return ( source.find(c) == 0 );
 333 }