BamAux.h

   1 // ***************************************************************************\r
   2 // BamAux.h (c) 2009 Derek Barnett, Michael Str�mberg\r
   3 // Marth Lab, Department of Biology, Boston College\r
   4 // All rights reserved.\r
   5 // ---------------------------------------------------------------------------\r
   6 // Last modified: 22 July 2010 (DB)\r
   7 // ---------------------------------------------------------------------------\r
   8 // Provides the basic constants, data structures, etc. for using BAM files\r
   9 // ***************************************************************************\r
  10 \r
  11 #ifndef BAMAUX_H\r
  12 #define BAMAUX_H\r
  13 \r
  14 // C inclues\r
  15 #include <cctype>\r
  16 #include <cstdio>\r
  17 #include <cstdlib>\r
  18 #include <cstring>\r
  19 \r
  20 // C++ includes\r
  21 #include <exception>\r
  22 #include <map>\r
  23 #include <string>\r
  24 #include <utility>\r
  25 #include <vector>\r
  26 \r
  27 #include <iostream>\r
  28 #include <typeinfo>\r
  29 \r
  30 // Platform-specific type definitions\r
  31 #ifndef BAMTOOLS_TYPES\r
  32 #define BAMTOOLS_TYPES\r
  33     #ifdef _MSC_VER\r
  34         typedef char                 int8_t;\r
  35         typedef unsigned char       uint8_t;\r
  36         typedef short               int16_t;\r
  37         typedef unsigned short     uint16_t;\r
  38         typedef int                 int32_t;\r
  39         typedef unsigned int       uint32_t;\r
  40         typedef long long           int64_t;\r
  41         typedef unsigned long long uint64_t;\r
  42     #else\r
  43         #include <stdint.h>\r
  44     #endif\r
  45 #endif // BAMTOOLS_TYPES\r
  46 \r
  47 namespace BamTools {\r
  48 \r
  49 // BAM constants\r
  50 const int BAM_CORE_SIZE   = 32;\r
  51 const int BAM_CMATCH      = 0;\r
  52 const int BAM_CINS        = 1;\r
  53 const int BAM_CDEL        = 2;\r
  54 const int BAM_CREF_SKIP   = 3;\r
  55 const int BAM_CSOFT_CLIP  = 4;\r
  56 const int BAM_CHARD_CLIP  = 5;\r
  57 const int BAM_CPAD        = 6;\r
  58 const int BAM_CIGAR_SHIFT = 4;\r
  59 const int BAM_CIGAR_MASK  = ((1 << BAM_CIGAR_SHIFT) - 1);\r
  60 \r
  61 // BAM index constants\r
  62 const int MAX_BIN           = 37450;    // =(8^6-1)/7+1\r
  63 const int BAM_MIN_CHUNK_GAP = 32768;\r
  64 const int BAM_LIDX_SHIFT    = 14;\r
  65 \r
  66 // Explicit variable sizes\r
  67 const int BT_SIZEOF_INT = 4;\r
  68 \r
  69 struct CigarOp;\r
  70 \r
  71 struct BamAlignment {\r
  72 \r
  73     // constructors & destructor\r
  74     public:\r
  75         BamAlignment(void);\r
  76         BamAlignment(const BamAlignment& other);\r
  77         ~BamAlignment(void);\r
  78 \r
  79     // Queries against alignment flags\r
  80     public:        \r
  81         bool IsDuplicate(void) const;           // Returns true if this read is a PCR duplicate       \r
  82         bool IsFailedQC(void) const;            // Returns true if this read failed quality control      \r
  83         bool IsFirstMate(void) const;           // Returns true if alignment is first mate on read        \r
  84         bool IsMapped(void) const;              // Returns true if alignment is mapped        \r
  85         bool IsMateMapped(void) const;          // Returns true if alignment's mate is mapped        \r
  86         bool IsMateReverseStrand(void) const;   // Returns true if alignment's mate mapped to reverse strand        \r
  87         bool IsPaired(void) const;              // Returns true if alignment part of paired-end read        \r
  88         bool IsPrimaryAlignment(void) const;    // Returns true if reported position is primary alignment       \r
  89         bool IsProperPair(void) const;          // Returns true if alignment is part of read that satisfied paired-end resolution     \r
  90         bool IsReverseStrand(void) const;       // Returns true if alignment mapped to reverse strand\r
  91         bool IsSecondMate(void) const;          // Returns true if alignment is second mate on read\r
  92 \r
  93     // Manipulate alignment flags\r
  94     public:        \r
  95         void SetIsDuplicate(bool ok);           // Sets "PCR duplicate" flag        \r
  96         void SetIsFailedQC(bool ok);            // Sets "failed quality control" flag        \r
  97         void SetIsFirstMate(bool ok);           // Sets "alignment is first mate" flag        \r
  98         void SetIsMateUnmapped(bool ok);        // Sets "alignment's mate is mapped" flag        \r
  99         void SetIsMateReverseStrand(bool ok);   // Sets "alignment's mate mapped to reverse strand" flag        \r
 100         void SetIsPaired(bool ok);              // Sets "alignment part of paired-end read" flag        \r
 101         void SetIsProperPair(bool ok);          // Sets "alignment is part of read that satisfied paired-end resolution" flag        \r
 102         void SetIsReverseStrand(bool ok);       // Sets "alignment mapped to reverse strand" flag        \r
 103         void SetIsSecondaryAlignment(bool ok);  // Sets "position is primary alignment" flag        \r
 104         void SetIsSecondMate(bool ok);          // Sets "alignment is second mate on read" flag        \r
 105         void SetIsUnmapped(bool ok);            // Sets "alignment is mapped" flag\r
 106 \r
 107     // Tag data access methods\r
 108     public:\r
 109         // generic tag data access methods \r
 110         bool GetTag(const std::string& tag, std::string& destination) const;    // access variable-length char or hex strings \r
 111         bool GetTag(const std::string& tag, uint32_t& destination) const;       // access unsigned integer data\r
 112         bool GetTag(const std::string& tag, int32_t& destination) const;        // access signed integer data\r
 113         bool GetTag(const std::string& tag, float& destination) const;          // access floating point data\r
 114         \r
 115         // specific tag data access methods - only remain here for legacy support\r
 116         bool GetEditDistance(uint8_t& editDistance) const;      // get "NM" tag data - contributed by Aaron Quinlan\r
 117         bool GetReadGroup(std::string& readGroup) const;        // get "RG" tag data\r
 118 \r
 119 \r
 120     // Additional data access methods\r
 121     public:\r
 122         int GetEndPosition(bool usePadded = false) const;       // calculates alignment end position, based on starting position and CIGAR operations\r
 123 \r
 124     // 'internal' utility methods \r
 125     private:\r
 126         static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);\r
 127 \r
 128     // Data members\r
 129     public:\r
 130         std::string  Name;              // Read name\r
 131         int32_t      Length;            // Query length\r
 132         std::string  QueryBases;        // 'Original' sequence (as reported from sequencing machine)\r
 133         std::string  AlignedBases;      // 'Aligned' sequence (includes any indels, padding, clipping)\r
 134         std::string  Qualities;         // FASTQ qualities (ASCII characters, not numeric values)\r
 135         std::string  TagData;           // Tag data (accessor methods will pull the requested information out)\r
 136         int32_t      RefID;             // ID number for reference sequence\r
 137         int32_t      Position;          // Position (0-based) where alignment starts\r
 138         uint16_t     Bin;               // Bin in BAM file where this alignment resides\r
 139         uint16_t     MapQuality;        // Mapping quality score\r
 140         uint32_t     AlignmentFlag;     // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate \r
 141         std::vector<CigarOp> CigarData; // CIGAR operations for this alignment\r
 142         int32_t      MateRefID;         // ID number for reference sequence where alignment's mate was aligned\r
 143         int32_t      MatePosition;      // Position (0-based) where alignment's mate starts\r
 144         int32_t      InsertSize;        // Mate-pair insert size\r
 145           \r
 146           \r
 147         struct BamAlignmentSupportData {\r
 148       \r
 149             // data members\r
 150             std::string AllCharData;\r
 151             uint32_t    BlockLength;\r
 152             uint32_t    NumCigarOperations;\r
 153             uint32_t    QueryNameLength;\r
 154             uint32_t    QuerySequenceLength;\r
 155             bool        HasCoreOnly;\r
 156             \r
 157             // constructor\r
 158             BamAlignmentSupportData(void)\r
 159                 : BlockLength(0)\r
 160                 , NumCigarOperations(0)\r
 161                 , QueryNameLength(0)\r
 162                 , QuerySequenceLength(0)\r
 163                 , HasCoreOnly(false)\r
 164             { }\r
 165         };\r
 166         \r
 167         BamAlignmentSupportData SupportData;  // Contains raw character data & lengths \r
 168 \r
 169     // Alignment flag query constants\r
 170     // Use the get/set methods above instead\r
 171     private:\r
 172         enum { PAIRED        = 1\r
 173              , PROPER_PAIR   = 2\r
 174              , UNMAPPED      = 4\r
 175              , MATE_UNMAPPED = 8\r
 176              , REVERSE       = 16\r
 177              , MATE_REVERSE  = 32\r
 178              , READ_1        = 64\r
 179              , READ_2        = 128\r
 180              , SECONDARY     = 256\r
 181              , QC_FAILED     = 512\r
 182              , DUPLICATE     = 1024 \r
 183              };\r
 184 };\r
 185 \r
 186 // ----------------------------------------------------------------\r
 187 // Auxiliary data structs & typedefs\r
 188 \r
 189 struct CigarOp {\r
 190   \r
 191     // data members\r
 192     char     Type;   // Operation type (MIDNSHP)\r
 193     uint32_t Length; // Operation length (number of bases)\r
 194     \r
 195     // constructor\r
 196     CigarOp(const char type = '\0', \r
 197             const uint32_t length = 0) \r
 198         : Type(type)\r
 199         , Length(length) \r
 200     { }\r
 201 };\r
 202 \r
 203 struct RefData {\r
 204    \r
 205     // data members\r
 206     std::string RefName;          // Name of reference sequence\r
 207     int32_t     RefLength;        // Length of reference sequence\r
 208     bool        RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence\r
 209     \r
 210     // constructor\r
 211     RefData(const int32_t& length = 0, \r
 212             bool ok = false)\r
 213         : RefLength(length)\r
 214         , RefHasAlignments(ok)\r
 215     { }\r
 216 };\r
 217 \r
 218 typedef std::vector<RefData>      RefVector;\r
 219 typedef std::vector<BamAlignment> BamAlignmentVector;\r
 220 \r
 221 struct BamRegion {\r
 222   \r
 223     // data members\r
 224     int LeftRefID;\r
 225     int LeftPosition;\r
 226     int RightRefID;\r
 227     int RightPosition;\r
 228     \r
 229     // constructor\r
 230     BamRegion(const int& leftID   = -1, \r
 231               const int& leftPos  = -1,\r
 232               const int& rightID  = -1,\r
 233               const int& rightPos = -1)\r
 234         : LeftRefID(leftID)\r
 235         , LeftPosition(leftPos)\r
 236         , RightRefID(rightID)\r
 237         , RightPosition(rightPos)\r
 238     { }\r
 239 };\r
 240 \r
 241 // ----------------------------------------------------------------\r
 242 // BamAlignment member methods\r
 243 \r
 244 // constructors & destructor\r
 245 inline \r
 246 BamAlignment::BamAlignment(void) { }\r
 247 \r
 248 inline \r
 249 BamAlignment::BamAlignment(const BamAlignment& other)\r
 250     : Name(other.Name)\r
 251     , Length(other.Length)\r
 252     , QueryBases(other.QueryBases)\r
 253     , AlignedBases(other.AlignedBases)\r
 254     , Qualities(other.Qualities)\r
 255     , TagData(other.TagData)\r
 256     , RefID(other.RefID)\r
 257     , Position(other.Position)\r
 258     , Bin(other.Bin)\r
 259     , MapQuality(other.MapQuality)\r
 260     , AlignmentFlag(other.AlignmentFlag)\r
 261     , CigarData(other.CigarData)\r
 262     , MateRefID(other.MateRefID)\r
 263     , MatePosition(other.MatePosition)\r
 264     , InsertSize(other.InsertSize)\r
 265     , SupportData(other.SupportData)\r
 266 { }\r
 267 \r
 268 inline \r
 269 BamAlignment::~BamAlignment(void) { }\r
 270 \r
 271 // Queries against alignment flags\r
 272 inline bool BamAlignment::IsDuplicate(void) const         { return ( (AlignmentFlag & DUPLICATE)     != 0 ); }\r
 273 inline bool BamAlignment::IsFailedQC(void) const          { return ( (AlignmentFlag & QC_FAILED)     != 0 ); }\r
 274 inline bool BamAlignment::IsFirstMate(void) const         { return ( (AlignmentFlag & READ_1)        != 0 ); }\r
 275 inline bool BamAlignment::IsMapped(void) const            { return ( (AlignmentFlag & UNMAPPED)      == 0 ); }\r
 276 inline bool BamAlignment::IsMateMapped(void) const        { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }\r
 277 inline bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }\r
 278 inline bool BamAlignment::IsPaired(void) const            { return ( (AlignmentFlag & PAIRED)        != 0 ); }\r
 279 inline bool BamAlignment::IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY)     == 0 ); }\r
 280 inline bool BamAlignment::IsProperPair(void) const        { return ( (AlignmentFlag & PROPER_PAIR)   != 0 ); }\r
 281 inline bool BamAlignment::IsReverseStrand(void) const     { return ( (AlignmentFlag & REVERSE)       != 0 ); }\r
 282 inline bool BamAlignment::IsSecondMate(void) const        { return ( (AlignmentFlag & READ_2)        != 0 ); }\r
 283 \r
 284 // Manipulate alignment flags \r
 285 inline void BamAlignment::SetIsDuplicate(bool ok)          { if (ok) AlignmentFlag |= DUPLICATE;     else AlignmentFlag &= ~DUPLICATE; }\r
 286 inline void BamAlignment::SetIsFailedQC(bool ok)           { if (ok) AlignmentFlag |= QC_FAILED;     else AlignmentFlag &= ~QC_FAILED; }\r
 287 inline void BamAlignment::SetIsFirstMate(bool ok)          { if (ok) AlignmentFlag |= READ_1;        else AlignmentFlag &= ~READ_1; }\r
 288 inline void BamAlignment::SetIsMateUnmapped(bool ok)       { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }\r
 289 inline void BamAlignment::SetIsMateReverseStrand(bool ok)  { if (ok) AlignmentFlag |= MATE_REVERSE;  else AlignmentFlag &= ~MATE_REVERSE; }\r
 290 inline void BamAlignment::SetIsPaired(bool ok)             { if (ok) AlignmentFlag |= PAIRED;        else AlignmentFlag &= ~PAIRED; }\r
 291 inline void BamAlignment::SetIsProperPair(bool ok)         { if (ok) AlignmentFlag |= PROPER_PAIR;   else AlignmentFlag &= ~PROPER_PAIR; }\r
 292 inline void BamAlignment::SetIsReverseStrand(bool ok)      { if (ok) AlignmentFlag |= REVERSE;       else AlignmentFlag &= ~REVERSE; }\r
 293 inline void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY;     else AlignmentFlag &= ~SECONDARY; }\r
 294 inline void BamAlignment::SetIsSecondMate(bool ok)         { if (ok) AlignmentFlag |= READ_2;        else AlignmentFlag &= ~READ_2; }\r
 295 inline void BamAlignment::SetIsUnmapped(bool ok)           { if (ok) AlignmentFlag |= UNMAPPED;      else AlignmentFlag &= ~UNMAPPED; }\r
 296 \r
 297 // calculates alignment end position, based on starting position and CIGAR operations\r
 298 inline \r
 299 int BamAlignment::GetEndPosition(bool usePadded) const {\r
 300 \r
 301     // initialize alignment end to starting position\r
 302     int alignEnd = Position;\r
 303 \r
 304     // iterate over cigar operations\r
 305     std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();\r
 306     std::vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();\r
 307     for ( ; cigarIter != cigarEnd; ++cigarIter) {\r
 308         const char cigarType = (*cigarIter).Type;\r
 309         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {\r
 310             alignEnd += (*cigarIter).Length;\r
 311         } \r
 312         else if ( usePadded && cigarType == 'I' ) {\r
 313             alignEnd += (*cigarIter).Length;\r
 314         }\r
 315     }\r
 316     return alignEnd;\r
 317 }\r
 318 \r
 319 // get "NM" tag data - contributed by Aaron Quinlan\r
 320 // stores data in 'editDistance', returns success/fail\r
 321 inline \r
 322 bool BamAlignment::GetEditDistance(uint8_t& editDistance) const {\r
 323 \r
 324     // make sure tag data exists\r
 325     if ( TagData.empty() ) return false;\r
 326 \r
 327     // localize the tag data\r
 328     char* pTagData = (char*)TagData.data();\r
 329     const unsigned int tagDataLen = TagData.size();\r
 330     unsigned int numBytesParsed = 0;\r
 331 \r
 332     bool foundEditDistanceTag = false;\r
 333     while ( numBytesParsed < tagDataLen ) {\r
 334 \r
 335         const char* pTagType = pTagData;\r
 336         const char* pTagStorageType = pTagData + 2;\r
 337         pTagData       += 3;\r
 338         numBytesParsed += 3;\r
 339 \r
 340         // check the current tag\r
 341         if ( strncmp(pTagType, "NM", 2) == 0 ) {\r
 342             foundEditDistanceTag = true;\r
 343             break;\r
 344         }\r
 345 \r
 346         // get the storage class and find the next tag\r
 347         if ( *pTagStorageType == '\0' ) return false;\r
 348         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 349         if ( *pTagData == '\0' ) return false;\r
 350     }\r
 351     // return if the edit distance tag was not present\r
 352     if ( !foundEditDistanceTag ) return false;\r
 353 \r
 354     // assign the editDistance value\r
 355     std::memcpy(&editDistance, pTagData, 1);\r
 356     return true;\r
 357 }\r
 358 \r
 359 // get "RG" tag data\r
 360 // stores data in 'readGroup', returns success/fail\r
 361 inline \r
 362 bool BamAlignment::GetReadGroup(std::string& readGroup) const {\r
 363 \r
 364     // make sure tag data exists\r
 365     if ( TagData.empty() ) return false;\r
 366 \r
 367     // localize the tag data\r
 368     char* pTagData = (char*)TagData.data();\r
 369     const unsigned int tagDataLen = TagData.size();\r
 370     unsigned int numBytesParsed = 0;\r
 371 \r
 372     bool foundReadGroupTag = false;\r
 373     while( numBytesParsed < tagDataLen ) {\r
 374 \r
 375         const char* pTagType = pTagData;\r
 376         const char* pTagStorageType = pTagData + 2;\r
 377         pTagData       += 3;\r
 378         numBytesParsed += 3;\r
 379 \r
 380         // check the current tag\r
 381         if ( std::strncmp(pTagType, "RG", 2) == 0 ) {\r
 382             foundReadGroupTag = true;\r
 383             break;\r
 384         }\r
 385 \r
 386         // get the storage class and find the next tag\r
 387         if ( *pTagStorageType == '\0' ) return false;\r
 388         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 389         if ( *pTagData == '\0' ) return false;\r
 390     }\r
 391 \r
 392     // return if the read group tag was not present\r
 393     if ( !foundReadGroupTag ) return false;\r
 394 \r
 395     // assign the read group\r
 396     const unsigned int readGroupLen = std::strlen(pTagData);\r
 397     readGroup.resize(readGroupLen);\r
 398     std::memcpy( (char*)readGroup.data(), pTagData, readGroupLen );\r
 399     return true;\r
 400 }\r
 401 \r
 402 inline\r
 403 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {\r
 404   \r
 405     // make sure tag data exists\r
 406     if ( TagData.empty() ) return false;\r
 407 \r
 408     // localize the tag data\r
 409     char* pTagData = (char*)TagData.data();\r
 410     const unsigned int tagDataLen = TagData.size();\r
 411     unsigned int numBytesParsed = 0;\r
 412 \r
 413     bool foundReadGroupTag = false;\r
 414     while ( numBytesParsed < tagDataLen ) {\r
 415 \r
 416         const char* pTagType = pTagData;\r
 417         const char* pTagStorageType = pTagData + 2;\r
 418         pTagData       += 3;\r
 419         numBytesParsed += 3;\r
 420 \r
 421         // check the current tag\r
 422         if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 ) {\r
 423             foundReadGroupTag = true;\r
 424             break;\r
 425         }\r
 426 \r
 427         // get the storage class and find the next tag\r
 428         if ( *pTagStorageType == '\0' ) return false; \r
 429         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 430         if ( *pTagData == '\0' ) return false;\r
 431     }\r
 432 \r
 433     // return if the read group tag was not present\r
 434     if ( !foundReadGroupTag ) return false;\r
 435 \r
 436     // assign the read group\r
 437     const unsigned int dataLen = std::strlen(pTagData);\r
 438     destination.resize(dataLen);\r
 439     std::memcpy( (char*)destination.data(), pTagData, dataLen );\r
 440     return true;\r
 441 }\r
 442 \r
 443 inline\r
 444 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {\r
 445   \r
 446     // make sure data exists\r
 447     if ( TagData.empty() ) return false;\r
 448 \r
 449     // clear out destination\r
 450     destination = 0;\r
 451 \r
 452     // localize the tag data\r
 453     char* pTagData = (char*)TagData.data();\r
 454     const unsigned int tagDataLen = TagData.size();\r
 455     unsigned int numBytesParsed = 0;\r
 456 \r
 457     int destinationLength = 0;    \r
 458     bool foundDesiredTag = false;\r
 459     while ( numBytesParsed < tagDataLen ) {\r
 460 \r
 461         const char* pTagType = pTagData;\r
 462         const char* pTagStorageType = pTagData + 2;\r
 463         pTagData       += 3;\r
 464         numBytesParsed += 3;\r
 465 \r
 466         // check the current tag\r
 467         if ( strncmp(pTagType, tag.c_str(), 2) == 0 ) {\r
 468             \r
 469             // determine actual length of data depending on tag type\r
 470             // this is necessary because some tags may be of variable byte-lengths (i.e. char or short)\r
 471             const char type = *pTagStorageType;\r
 472             switch(type) {\r
 473 \r
 474                 // 1 byte data\r
 475                 case 'A':\r
 476                 case 'c':\r
 477                 case 'C':\r
 478                     destinationLength = 1;\r
 479                     break;\r
 480 \r
 481                 // 2 byte data\r
 482                 case 's':\r
 483                 case 'S':\r
 484                     destinationLength = 2;\r
 485                     break;\r
 486 \r
 487                 // 4 byte data\r
 488                 case 'i':\r
 489                 case 'I':\r
 490                     destinationLength = 4;\r
 491                     break;\r
 492 \r
 493                 // unsupported type for integer destination (float & var-length strings)\r
 494                 case 'f':\r
 495                 case 'Z':\r
 496                 case 'H':\r
 497                     printf("ERROR: Cannot store tag of type %c in integer destination\n", type);\r
 498                     return false;\r
 499 \r
 500                 // unknown tag type\r
 501                 default:\r
 502                     printf("ERROR: Unknown tag storage class encountered: [%c]\n", *pTagData);\r
 503                     return false;\r
 504             }\r
 505             \r
 506             foundDesiredTag = true;\r
 507             break;\r
 508         }\r
 509 \r
 510         // get the storage class and find the next tag\r
 511         if ( *pTagStorageType == '\0' ) return false;\r
 512         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 513         if ( *pTagData == '\0' ) return false;\r
 514     }\r
 515     // return if the edit distance tag was not present\r
 516     if ( !foundDesiredTag ) return false; \r
 517 \r
 518     // assign the editDistance value\r
 519     std::memcpy(&destination, pTagData, destinationLength);\r
 520     return true;\r
 521 }\r
 522 \r
 523 inline\r
 524 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {\r
 525     return GetTag(tag, (uint32_t&)destination);\r
 526 }\r
 527 \r
 528 inline\r
 529 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {\r
 530   \r
 531     // make sure data exists\r
 532     if ( TagData.empty() ) return false;\r
 533 \r
 534     // clear out destination\r
 535     destination = 0.0;\r
 536 \r
 537     // localize the tag data\r
 538     char* pTagData = (char*)TagData.data();\r
 539     const unsigned int tagDataLen = TagData.size();\r
 540     unsigned int numBytesParsed = 0;\r
 541 \r
 542     int destinationLength = 0;\r
 543     bool foundDesiredTag = false;\r
 544     while( numBytesParsed < tagDataLen ) {\r
 545 \r
 546         const char* pTagType = pTagData;\r
 547         const char* pTagStorageType = pTagData + 2;\r
 548         pTagData       += 3;\r
 549         numBytesParsed += 3;\r
 550 \r
 551         // check the current tag\r
 552         if ( strncmp(pTagType, tag.c_str(), 2) == 0 ) {\r
 553           \r
 554             // determine actual length of data depending on tag type\r
 555             // this is necessary because some tags may be of variable byte-lengths (i.e. char or short)\r
 556             const char type = *pTagStorageType;\r
 557             switch(type) {\r
 558 \r
 559                 // 1 byte data\r
 560                 case 'A':\r
 561                 case 'c':\r
 562                 case 'C':\r
 563                     destinationLength = 1;\r
 564                     break;\r
 565 \r
 566                 // 2 byte data\r
 567                 case 's':\r
 568                 case 'S':\r
 569                     destinationLength = 2;\r
 570                     break;\r
 571 \r
 572                 // 4 byte data\r
 573                 case 'f':\r
 574                 case 'i':\r
 575                 case 'I':\r
 576                     destinationLength = 4;\r
 577                     break;\r
 578                 \r
 579                 // unsupported type (var-length strings)\r
 580                 case 'Z':\r
 581                 case 'H':\r
 582                     printf("ERROR: Cannot store tag of type %c in integer destination\n", type);\r
 583                     return false;\r
 584 \r
 585                 // unknown tag type\r
 586                 default:\r
 587                     printf("ERROR: Unknown tag storage class encountered: [%c]\n", *pTagData);\r
 588                     return false;\r
 589             }\r
 590             \r
 591             foundDesiredTag = true;\r
 592             break;\r
 593         }\r
 594 \r
 595         // get the storage class and find the next tag\r
 596         if ( *pTagStorageType == '\0' ) return false;\r
 597         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 598         if ( *pTagData == '\0' ) return false;\r
 599     }\r
 600     // return if the edit distance tag was not present\r
 601     if ( !foundDesiredTag ) return false; \r
 602 \r
 603     // assign the editDistance value\r
 604     std::memcpy(&destination, pTagData, destinationLength);\r
 605     return true;\r
 606 }\r
 607 \r
 608 inline\r
 609 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {\r
 610     \r
 611     switch(storageType) {\r
 612 \r
 613         case 'A':\r
 614         case 'c':\r
 615         case 'C':\r
 616             ++numBytesParsed;\r
 617             ++pTagData;\r
 618             break;\r
 619 \r
 620         case 's':\r
 621         case 'S':\r
 622             numBytesParsed += 2;\r
 623             pTagData       += 2;\r
 624             break;\r
 625 \r
 626         case 'f':\r
 627         case 'i':\r
 628         case 'I':\r
 629             numBytesParsed += 4;\r
 630             pTagData       += 4;\r
 631             break;\r
 632 \r
 633         case 'Z':\r
 634         case 'H':\r
 635             while(*pTagData) {\r
 636                 ++numBytesParsed;\r
 637                 ++pTagData;\r
 638             }\r
 639         // ---------------------------\r
 640         // Added: 3-25-2010 DWB\r
 641         // Contributed: ARQ\r
 642         // Fixed: error parsing variable length tag data\r
 643             ++pTagData;\r
 644         // ---------------------------\r
 645             break;\r
 646 \r
 647         default: \r
 648             // error case\r
 649             printf("ERROR: Unknown tag storage class encountered: [%c]\n", *pTagData);\r
 650             return false;\r
 651     }\r
 652     \r
 653     // return success\r
 654     return true;\r
 655 }\r
 656 \r
 657 // ----------------------------------------------------------------\r
 658 // Added: 3-35-2010 DWB\r
 659 // Fixed: Routines to provide endian-correctness\r
 660 // ----------------------------------------------------------------\r
 661 \r
 662 // returns true if system is big endian\r
 663 inline bool SystemIsBigEndian(void) {\r
 664    const uint16_t one = 0x0001;\r
 665    return ((*(char*) &one) == 0 );\r
 666 }\r
 667 \r
 668 // swaps endianness of 16-bit value 'in place'\r
 669 inline void SwapEndian_16(int16_t& x) {\r
 670     x = ((x >> 8) | (x << 8));\r
 671 }\r
 672 \r
 673 inline void SwapEndian_16(uint16_t& x) {\r
 674     x = ((x >> 8) | (x << 8));\r
 675 }\r
 676 \r
 677 // swaps endianness of 32-bit value 'in-place'\r
 678 inline void SwapEndian_32(int32_t& x) {\r
 679     x = ( (x >> 24) | \r
 680          ((x << 8) & 0x00FF0000) | \r
 681          ((x >> 8) & 0x0000FF00) | \r
 682           (x << 24)\r
 683         );\r
 684 }\r
 685 \r
 686 inline void SwapEndian_32(uint32_t& x) {\r
 687     x = ( (x >> 24) | \r
 688          ((x << 8) & 0x00FF0000) | \r
 689          ((x >> 8) & 0x0000FF00) | \r
 690           (x << 24)\r
 691         );\r
 692 }\r
 693 \r
 694 // swaps endianness of 64-bit value 'in-place'\r
 695 inline void SwapEndian_64(int64_t& x) {\r
 696     x = ( (x >> 56) | \r
 697          ((x << 40) & 0x00FF000000000000ll) |\r
 698          ((x << 24) & 0x0000FF0000000000ll) |\r
 699          ((x << 8)  & 0x000000FF00000000ll) |\r
 700          ((x >> 8)  & 0x00000000FF000000ll) |\r
 701          ((x >> 24) & 0x0000000000FF0000ll) |\r
 702          ((x >> 40) & 0x000000000000FF00ll) |\r
 703           (x << 56)\r
 704         );\r
 705 }\r
 706 \r
 707 inline void SwapEndian_64(uint64_t& x) {\r
 708     x = ( (x >> 56) | \r
 709          ((x << 40) & 0x00FF000000000000ll) |\r
 710          ((x << 24) & 0x0000FF0000000000ll) |\r
 711          ((x << 8)  & 0x000000FF00000000ll) |\r
 712          ((x >> 8)  & 0x00000000FF000000ll) |\r
 713          ((x >> 24) & 0x0000000000FF0000ll) |\r
 714          ((x >> 40) & 0x000000000000FF00ll) |\r
 715           (x << 56)\r
 716         );\r
 717 }\r
 718 \r
 719 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)\r
 720 inline void SwapEndian_16p(char* data) {\r
 721     uint16_t& value = (uint16_t&)*data; \r
 722     SwapEndian_16(value);\r
 723 }\r
 724 \r
 725 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)\r
 726 inline void SwapEndian_32p(char* data) {\r
 727     uint32_t& value = (uint32_t&)*data; \r
 728     SwapEndian_32(value);\r
 729 }\r
 730 \r
 731 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)\r
 732 inline void SwapEndian_64p(char* data) {\r
 733     uint64_t& value = (uint64_t&)*data; \r
 734     SwapEndian_64(value);\r
 735 }\r
 736 \r
 737 } // namespace BamTools\r
 738 \r
 739 #endif // BAMAUX_H\r