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