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