src/api/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: 16 September 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 <fstream>\r
  23 #include <iostream>\r
  24 #include <map>\r
  25 #include <string>\r
  26 #include <utility>\r
  27 #include <vector>\r
  28 \r
  29 // Platform-specific type definitions\r
  30 #ifndef BAMTOOLS_TYPES\r
  31 #define BAMTOOLS_TYPES\r
  32     #ifdef _MSC_VER\r
  33         typedef char                 int8_t;\r
  34         typedef unsigned char       uint8_t;\r
  35         typedef short               int16_t;\r
  36         typedef unsigned short     uint16_t;\r
  37         typedef int                 int32_t;\r
  38         typedef unsigned int       uint32_t;\r
  39         typedef long long           int64_t;\r
  40         typedef unsigned long long uint64_t;\r
  41     #else\r
  42         #include <stdint.h>\r
  43     #endif\r
  44 #endif // BAMTOOLS_TYPES\r
  45 \r
  46 namespace BamTools {\r
  47 \r
  48 // BAM constants\r
  49 const int BAM_CIGAR_SHIFT = 4;\r
  50 const int BAM_CIGAR_MASK  = ((1 << BAM_CIGAR_SHIFT) - 1);\r
  51 const int BAM_CORE_SIZE   = 32;\r
  52 \r
  53 struct CigarOp;\r
  54 \r
  55 struct BamAlignment {\r
  56 \r
  57     // constructors & destructor\r
  58     public:\r
  59         BamAlignment(void);\r
  60         BamAlignment(const BamAlignment& other);\r
  61         ~BamAlignment(void);\r
  62 \r
  63     // Queries against alignment flags\r
  64     public:        \r
  65         bool IsDuplicate(void) const;           // Returns true if this read is a PCR duplicate       \r
  66         bool IsFailedQC(void) const;            // Returns true if this read failed quality control      \r
  67         bool IsFirstMate(void) const;           // Returns true if alignment is first mate on read        \r
  68         bool IsMapped(void) const;              // Returns true if alignment is mapped        \r
  69         bool IsMateMapped(void) const;          // Returns true if alignment's mate is mapped        \r
  70         bool IsMateReverseStrand(void) const;   // Returns true if alignment's mate mapped to reverse strand        \r
  71         bool IsPaired(void) const;              // Returns true if alignment part of paired-end read        \r
  72         bool IsPrimaryAlignment(void) const;    // Returns true if reported position is primary alignment       \r
  73         bool IsProperPair(void) const;          // Returns true if alignment is part of read that satisfied paired-end resolution     \r
  74         bool IsReverseStrand(void) const;       // Returns true if alignment mapped to reverse strand\r
  75         bool IsSecondMate(void) const;          // Returns true if alignment is second mate on read\r
  76 \r
  77     // Manipulate alignment flags\r
  78     public:        \r
  79         void SetIsDuplicate(bool ok);           // Sets "PCR duplicate" flag        \r
  80         void SetIsFailedQC(bool ok);            // Sets "failed quality control" flag        \r
  81         void SetIsFirstMate(bool ok);           // Sets "alignment is first mate" flag        \r
  82         void SetIsMateUnmapped(bool ok);        // Sets "alignment's mate is mapped" flag        \r
  83         void SetIsMateReverseStrand(bool ok);   // Sets "alignment's mate mapped to reverse strand" flag        \r
  84         void SetIsPaired(bool ok);              // Sets "alignment part of paired-end read" flag        \r
  85         void SetIsProperPair(bool ok);          // Sets "alignment is part of read that satisfied paired-end resolution" flag        \r
  86         void SetIsReverseStrand(bool ok);       // Sets "alignment mapped to reverse strand" flag        \r
  87         void SetIsSecondaryAlignment(bool ok);  // Sets "position is primary alignment" flag        \r
  88         void SetIsSecondMate(bool ok);          // Sets "alignment is second mate on read" flag        \r
  89         void SetIsUnmapped(bool ok);            // Sets "alignment is mapped" flag\r
  90 \r
  91     // Tag data access methods\r
  92     public:\r
  93         // -------------------------------------------------------------------------------------\r
  94         // N.B. - The following tag-modifying methods may not be used on BamAlignments fetched\r
  95         // using BamReader::GetNextAlignmentCore().  Attempting to use them will not result in \r
  96         // error message (to keep output clean) but will ALWAYS return false.  Only user-\r
  97         // generated BamAlignments or those retrieved using BamReader::GetNextAlignment() are valid.\r
  98 \r
  99         // add tag data (create new TAG entry with TYPE and VALUE)\r
 100         // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details\r
 101         // returns true if new data added, false if error or TAG already exists\r
 102         // N.B. - will NOT modify existing tag. Use EditTag() instead\r
 103         bool AddTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H\r
 104         bool AddTag(const std::string& tag, const std::string& type, const uint32_t& value);    // type must be A or i\r
 105         bool AddTag(const std::string& tag, const std::string& type, const int32_t& value);     // type must be A or i\r
 106         bool AddTag(const std::string& tag, const std::string& type, const float& value);       // type must be A, i, or f\r
 107         \r
 108         // edit tag data (sets existing TAG with TYPE to VALUE or adds new TAG if not already present)\r
 109         // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details\r
 110         // returns true if edit was successfaul, false if error\r
 111         bool EditTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H\r
 112         bool EditTag(const std::string& tag, const std::string& type, const uint32_t& value);    // type must be A or i\r
 113         bool EditTag(const std::string& tag, const std::string& type, const int32_t& value);     // type must be A or i\r
 114         bool EditTag(const std::string& tag, const std::string& type, const float& value);       // type must be A, i, or f\r
 115 \r
 116         // specific tag data access methods - these only remain for legacy support\r
 117         bool GetEditDistance(uint32_t& editDistance) const; // get "NM" tag data (implemented as GetTag("NM", editDistance))\r
 118         bool GetReadGroup(std::string& readGroup) const;    // get "RG" tag data (implemented as GetTag("RG", readGroup)) \r
 119         \r
 120         // generic tag data access methods \r
 121         bool GetTag(const std::string& tag, std::string& destination) const;    // access variable-length char or hex strings \r
 122         bool GetTag(const std::string& tag, uint32_t& destination) const;       // access unsigned integer data\r
 123         bool GetTag(const std::string& tag, int32_t& destination) const;        // access signed integer data\r
 124         bool GetTag(const std::string& tag, float& destination) const;          // access floating point data\r
 125         \r
 126         // remove tag data\r
 127         // returns true if removal was successful, false if error\r
 128         // N.B. - returns false if TAG does not exist (no removal can occur)\r
 129         bool RemoveTag(const std::string& tag);\r
 130 \r
 131     // Additional data access methods\r
 132     public:\r
 133         // calculates alignment end position, based on starting position and CIGAR operations\r
 134         // @zeroBased - if true, returns 0-based coordinate; else returns 1-based\r
 135         int GetEndPosition(bool usePadded = false, bool zeroBased = true) const;  \r
 136 \r
 137     // 'internal' utility methods \r
 138     private:\r
 139         static bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed);\r
 140         static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);\r
 141 \r
 142     // Data members\r
 143     public:\r
 144         std::string Name;              // Read name\r
 145         int32_t     Length;            // Query length\r
 146         std::string QueryBases;        // 'Original' sequence (as reported from sequencing machine)\r
 147         std::string AlignedBases;      // 'Aligned' sequence (includes any indels, padding, clipping)\r
 148         std::string Qualities;         // FASTQ qualities (ASCII characters, not numeric values)\r
 149         std::string TagData;           // Tag data (accessor methods will pull the requested information out)\r
 150         int32_t     RefID;             // ID number for reference sequence\r
 151         int32_t     Position;          // Position (0-based) where alignment starts\r
 152         uint16_t    Bin;               // Bin in BAM file where this alignment resides\r
 153         uint16_t    MapQuality;        // Mapping quality score\r
 154         uint32_t    AlignmentFlag;     // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate \r
 155         std::vector<CigarOp> CigarData; // CIGAR operations for this alignment\r
 156         int32_t     MateRefID;         // ID number for reference sequence where alignment's mate was aligned\r
 157         int32_t     MatePosition;      // Position (0-based) where alignment's mate starts\r
 158         int32_t     InsertSize;        // Mate-pair insert size\r
 159           \r
 160     // internal data\r
 161     private:\r
 162         struct BamAlignmentSupportData {\r
 163       \r
 164             // data members\r
 165             std::string AllCharData;\r
 166             uint32_t    BlockLength;\r
 167             uint32_t    NumCigarOperations;\r
 168             uint32_t    QueryNameLength;\r
 169             uint32_t    QuerySequenceLength;\r
 170             bool        HasCoreOnly;\r
 171             \r
 172             // constructor\r
 173             BamAlignmentSupportData(void)\r
 174                 : BlockLength(0)\r
 175                 , NumCigarOperations(0)\r
 176                 , QueryNameLength(0)\r
 177                 , QuerySequenceLength(0)\r
 178                 , HasCoreOnly(false)\r
 179             { }\r
 180         };\r
 181         \r
 182         // contains raw character data & lengths\r
 183         BamAlignmentSupportData SupportData;   \r
 184         \r
 185         // allow these classes access to BamAlignment private members (SupportData)\r
 186         // but client code should not need to touch this data\r
 187         friend class BamReader;\r
 188         friend class BamWriter;\r
 189 \r
 190     // Alignment flag query constants\r
 191     // Use the get/set methods above instead\r
 192     private:\r
 193         enum { PAIRED        = 1\r
 194              , PROPER_PAIR   = 2\r
 195              , UNMAPPED      = 4\r
 196              , MATE_UNMAPPED = 8\r
 197              , REVERSE       = 16\r
 198              , MATE_REVERSE  = 32\r
 199              , READ_1        = 64\r
 200              , READ_2        = 128\r
 201              , SECONDARY     = 256\r
 202              , QC_FAILED     = 512\r
 203              , DUPLICATE     = 1024 \r
 204              };\r
 205 };\r
 206 \r
 207 // ----------------------------------------------------------------\r
 208 // Auxiliary data structs & typedefs\r
 209 \r
 210 struct CigarOp {\r
 211   \r
 212     // data members\r
 213     char     Type;   // Operation type (MIDNSHP)\r
 214     uint32_t Length; // Operation length (number of bases)\r
 215     \r
 216     // constructor\r
 217     CigarOp(const char type = '\0', \r
 218             const uint32_t length = 0) \r
 219         : Type(type)\r
 220         , Length(length) \r
 221     { }\r
 222 };\r
 223 \r
 224 struct RefData {\r
 225    \r
 226     // data members\r
 227     std::string RefName;          // Name of reference sequence\r
 228     int32_t     RefLength;        // Length of reference sequence\r
 229     bool        RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence\r
 230     \r
 231     // constructor\r
 232     RefData(const int32_t& length = 0, \r
 233             bool ok = false)\r
 234         : RefLength(length)\r
 235         , RefHasAlignments(ok)\r
 236     { }\r
 237 };\r
 238 \r
 239 typedef std::vector<RefData>      RefVector;\r
 240 typedef std::vector<BamAlignment> BamAlignmentVector;\r
 241 \r
 242 struct BamRegion {\r
 243   \r
 244     // data members\r
 245     int LeftRefID;\r
 246     int LeftPosition;\r
 247     int RightRefID;\r
 248     int RightPosition;\r
 249     \r
 250     // constructor\r
 251     BamRegion(const int& leftID   = -1, \r
 252               const int& leftPos  = -1,\r
 253               const int& rightID  = -1,\r
 254               const int& rightPos = -1)\r
 255         : LeftRefID(leftID)\r
 256         , LeftPosition(leftPos)\r
 257         , RightRefID(rightID)\r
 258         , RightPosition(rightPos)\r
 259     { }\r
 260     \r
 261     // member functions\r
 262     void clear(void) { LeftRefID = -1; LeftPosition = -1; RightRefID = -1; RightPosition = -1; }\r
 263     bool isLeftBoundSpecified(void) const { return ( LeftRefID != -1 && LeftPosition != -1 ); }\r
 264     bool isNull(void) const { return ( !isLeftBoundSpecified() && !isRightBoundSpecified() ); }\r
 265     bool isRightBoundSpecified(void) const { return ( RightRefID != -1 && RightPosition != -1 ); }\r
 266 };\r
 267 \r
 268 // ----------------------------------------------------------------\r
 269 // Added: 3-35-2010 DWB\r
 270 // Fixed: Routines to provide endian-correctness\r
 271 // ----------------------------------------------------------------\r
 272 \r
 273 // returns true if system is big endian\r
 274 inline bool SystemIsBigEndian(void) {\r
 275    const uint16_t one = 0x0001;\r
 276    return ((*(char*) &one) == 0 );\r
 277 }\r
 278 \r
 279 // swaps endianness of 16-bit value 'in place'\r
 280 inline void SwapEndian_16(int16_t& x) {\r
 281     x = ((x >> 8) | (x << 8));\r
 282 }\r
 283 \r
 284 inline void SwapEndian_16(uint16_t& x) {\r
 285     x = ((x >> 8) | (x << 8));\r
 286 }\r
 287 \r
 288 // swaps endianness of 32-bit value 'in-place'\r
 289 inline void SwapEndian_32(int32_t& x) {\r
 290     x = ( (x >> 24) | \r
 291          ((x << 8) & 0x00FF0000) | \r
 292          ((x >> 8) & 0x0000FF00) | \r
 293           (x << 24)\r
 294         );\r
 295 }\r
 296 \r
 297 inline void SwapEndian_32(uint32_t& x) {\r
 298     x = ( (x >> 24) | \r
 299          ((x << 8) & 0x00FF0000) | \r
 300          ((x >> 8) & 0x0000FF00) | \r
 301           (x << 24)\r
 302         );\r
 303 }\r
 304 \r
 305 // swaps endianness of 64-bit value 'in-place'\r
 306 inline void SwapEndian_64(int64_t& x) {\r
 307     x = ( (x >> 56) | \r
 308          ((x << 40) & 0x00FF000000000000ll) |\r
 309          ((x << 24) & 0x0000FF0000000000ll) |\r
 310          ((x << 8)  & 0x000000FF00000000ll) |\r
 311          ((x >> 8)  & 0x00000000FF000000ll) |\r
 312          ((x >> 24) & 0x0000000000FF0000ll) |\r
 313          ((x >> 40) & 0x000000000000FF00ll) |\r
 314           (x << 56)\r
 315         );\r
 316 }\r
 317 \r
 318 inline void SwapEndian_64(uint64_t& x) {\r
 319     x = ( (x >> 56) | \r
 320          ((x << 40) & 0x00FF000000000000ll) |\r
 321          ((x << 24) & 0x0000FF0000000000ll) |\r
 322          ((x << 8)  & 0x000000FF00000000ll) |\r
 323          ((x >> 8)  & 0x00000000FF000000ll) |\r
 324          ((x >> 24) & 0x0000000000FF0000ll) |\r
 325          ((x >> 40) & 0x000000000000FF00ll) |\r
 326           (x << 56)\r
 327         );\r
 328 }\r
 329 \r
 330 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)\r
 331 inline void SwapEndian_16p(char* data) {\r
 332     uint16_t& value = (uint16_t&)*data; \r
 333     SwapEndian_16(value);\r
 334 }\r
 335 \r
 336 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)\r
 337 inline void SwapEndian_32p(char* data) {\r
 338     uint32_t& value = (uint32_t&)*data; \r
 339     SwapEndian_32(value);\r
 340 }\r
 341 \r
 342 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)\r
 343 inline void SwapEndian_64p(char* data) {\r
 344     uint64_t& value = (uint64_t&)*data; \r
 345     SwapEndian_64(value);\r
 346 }\r
 347 \r
 348 inline bool FileExists(const std::string& filename) {\r
 349     std::ifstream f(filename.c_str(), std::ifstream::in);\r
 350     return !f.fail();\r
 351 }\r
 352 \r
 353 // ----------------------------------------------------------------\r
 354 // BamAlignment member methods\r
 355 \r
 356 // constructors & destructor\r
 357 inline BamAlignment::BamAlignment(void) { }\r
 358 \r
 359 inline BamAlignment::BamAlignment(const BamAlignment& other)\r
 360     : Name(other.Name)\r
 361     , Length(other.Length)\r
 362     , QueryBases(other.QueryBases)\r
 363     , AlignedBases(other.AlignedBases)\r
 364     , Qualities(other.Qualities)\r
 365     , TagData(other.TagData)\r
 366     , RefID(other.RefID)\r
 367     , Position(other.Position)\r
 368     , Bin(other.Bin)\r
 369     , MapQuality(other.MapQuality)\r
 370     , AlignmentFlag(other.AlignmentFlag)\r
 371     , CigarData(other.CigarData)\r
 372     , MateRefID(other.MateRefID)\r
 373     , MatePosition(other.MatePosition)\r
 374     , InsertSize(other.InsertSize)\r
 375     , SupportData(other.SupportData)\r
 376 { }\r
 377 \r
 378 inline BamAlignment::~BamAlignment(void) { }\r
 379 \r
 380 // Queries against alignment flags\r
 381 inline bool BamAlignment::IsDuplicate(void) const         { return ( (AlignmentFlag & DUPLICATE)     != 0 ); }\r
 382 inline bool BamAlignment::IsFailedQC(void) const          { return ( (AlignmentFlag & QC_FAILED)     != 0 ); }\r
 383 inline bool BamAlignment::IsFirstMate(void) const         { return ( (AlignmentFlag & READ_1)        != 0 ); }\r
 384 inline bool BamAlignment::IsMapped(void) const            { return ( (AlignmentFlag & UNMAPPED)      == 0 ); }\r
 385 inline bool BamAlignment::IsMateMapped(void) const        { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }\r
 386 inline bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }\r
 387 inline bool BamAlignment::IsPaired(void) const            { return ( (AlignmentFlag & PAIRED)        != 0 ); }\r
 388 inline bool BamAlignment::IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY)     == 0 ); }\r
 389 inline bool BamAlignment::IsProperPair(void) const        { return ( (AlignmentFlag & PROPER_PAIR)   != 0 ); }\r
 390 inline bool BamAlignment::IsReverseStrand(void) const     { return ( (AlignmentFlag & REVERSE)       != 0 ); }\r
 391 inline bool BamAlignment::IsSecondMate(void) const        { return ( (AlignmentFlag & READ_2)        != 0 ); }\r
 392 \r
 393 // Manipulate alignment flags \r
 394 inline void BamAlignment::SetIsDuplicate(bool ok)          { if (ok) AlignmentFlag |= DUPLICATE;     else AlignmentFlag &= ~DUPLICATE; }\r
 395 inline void BamAlignment::SetIsFailedQC(bool ok)           { if (ok) AlignmentFlag |= QC_FAILED;     else AlignmentFlag &= ~QC_FAILED; }\r
 396 inline void BamAlignment::SetIsFirstMate(bool ok)          { if (ok) AlignmentFlag |= READ_1;        else AlignmentFlag &= ~READ_1; }\r
 397 inline void BamAlignment::SetIsMateUnmapped(bool ok)       { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }\r
 398 inline void BamAlignment::SetIsMateReverseStrand(bool ok)  { if (ok) AlignmentFlag |= MATE_REVERSE;  else AlignmentFlag &= ~MATE_REVERSE; }\r
 399 inline void BamAlignment::SetIsPaired(bool ok)             { if (ok) AlignmentFlag |= PAIRED;        else AlignmentFlag &= ~PAIRED; }\r
 400 inline void BamAlignment::SetIsProperPair(bool ok)         { if (ok) AlignmentFlag |= PROPER_PAIR;   else AlignmentFlag &= ~PROPER_PAIR; }\r
 401 inline void BamAlignment::SetIsReverseStrand(bool ok)      { if (ok) AlignmentFlag |= REVERSE;       else AlignmentFlag &= ~REVERSE; }\r
 402 inline void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY;     else AlignmentFlag &= ~SECONDARY; }\r
 403 inline void BamAlignment::SetIsSecondMate(bool ok)         { if (ok) AlignmentFlag |= READ_2;        else AlignmentFlag &= ~READ_2; }\r
 404 inline void BamAlignment::SetIsUnmapped(bool ok)           { if (ok) AlignmentFlag |= UNMAPPED;      else AlignmentFlag &= ~UNMAPPED; }\r
 405 \r
 406 // calculates alignment end position, based on starting position and CIGAR operations\r
 407 inline \r
 408 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {\r
 409 \r
 410     // initialize alignment end to starting position\r
 411     int alignEnd = Position;\r
 412 \r
 413     // iterate over cigar operations\r
 414     std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();\r
 415     std::vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();\r
 416     for ( ; cigarIter != cigarEnd; ++cigarIter) {\r
 417         const char cigarType = (*cigarIter).Type;\r
 418         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {\r
 419             alignEnd += (*cigarIter).Length;\r
 420         } \r
 421         else if ( usePadded && cigarType == 'I' ) {\r
 422             alignEnd += (*cigarIter).Length;\r
 423         }\r
 424     }\r
 425     \r
 426     // adjust for zeroBased, if necessary\r
 427     if (zeroBased) \r
 428         return alignEnd - 1;\r
 429     else \r
 430         return alignEnd;\r
 431 }\r
 432 \r
 433 inline\r
 434 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {\r
 435   \r
 436     if ( SupportData.HasCoreOnly ) return false;\r
 437     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 438     if ( type != "Z" && type != "H" ) return false;\r
 439   \r
 440     // localize the tag data\r
 441     char* pTagData = (char*)TagData.data();\r
 442     const unsigned int tagDataLength = TagData.size();\r
 443     unsigned int numBytesParsed = 0;\r
 444     \r
 445     // if tag already exists, return false\r
 446     // use EditTag explicitly instead\r
 447     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;\r
 448   \r
 449     // otherwise, copy tag data to temp buffer\r
 450     std::string newTag = tag + type + value;\r
 451     const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term\r
 452     char originalTagData[newTagDataLength];\r
 453     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term\r
 454     \r
 455     // append newTag\r
 456     strcat(originalTagData + tagDataLength, newTag.data());  // removes original null-term, appends newTag + null-term\r
 457     \r
 458     // store temp buffer back in TagData\r
 459     const char* newTagData = (const char*)originalTagData;\r
 460     TagData.assign(newTagData, newTagDataLength);\r
 461     \r
 462     // return success\r
 463     return true;\r
 464 }\r
 465 \r
 466 inline\r
 467 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {\r
 468   \r
 469     if ( SupportData.HasCoreOnly ) return false;\r
 470     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 471     if ( type == "f" || type == "Z" || type == "H" ) return false;\r
 472   \r
 473     // localize the tag data\r
 474     char* pTagData = (char*)TagData.data();\r
 475     const unsigned int tagDataLength = TagData.size();\r
 476     unsigned int numBytesParsed = 0;\r
 477     \r
 478     // if tag already exists, return false\r
 479     // use EditTag explicitly instead\r
 480     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;\r
 481   \r
 482     // otherwise, convert value to string\r
 483     union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;\r
 484     un.value = value;\r
 485 \r
 486     // copy original tag data to temp buffer\r
 487     std::string newTag = tag + type;\r
 488     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer\r
 489     char originalTagData[newTagDataLength];\r
 490     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term\r
 491     \r
 492     // append newTag\r
 493     strcat(originalTagData + tagDataLength, newTag.data());\r
 494     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));\r
 495     \r
 496     // store temp buffer back in TagData\r
 497     const char* newTagData = (const char*)originalTagData;\r
 498     TagData.assign(newTagData, newTagDataLength);\r
 499     \r
 500     // return success\r
 501     return true;\r
 502 }\r
 503 \r
 504 inline\r
 505 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {\r
 506     return AddTag(tag, type, (const uint32_t&)value);\r
 507 }\r
 508 \r
 509 inline\r
 510 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {\r
 511   \r
 512     if ( SupportData.HasCoreOnly ) return false;\r
 513     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 514     if ( type == "Z" || type == "H" ) return false;\r
 515   \r
 516     // localize the tag data\r
 517     char* pTagData = (char*)TagData.data();\r
 518     const unsigned int tagDataLength = TagData.size();\r
 519     unsigned int numBytesParsed = 0;\r
 520     \r
 521     // if tag already exists, return false\r
 522     // use EditTag explicitly instead\r
 523     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;\r
 524   \r
 525     // otherwise, convert value to string\r
 526     union { float value; char valueBuffer[sizeof(float)]; } un;\r
 527     un.value = value;\r
 528 \r
 529     // copy original tag data to temp buffer\r
 530     std::string newTag = tag + type;\r
 531     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float\r
 532     char originalTagData[newTagDataLength];\r
 533     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term\r
 534     \r
 535     // append newTag\r
 536     strcat(originalTagData + tagDataLength, newTag.data());\r
 537     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));\r
 538     \r
 539     // store temp buffer back in TagData\r
 540     const char* newTagData = (const char*)originalTagData;\r
 541     TagData.assign(newTagData, newTagDataLength);\r
 542     \r
 543     // return success\r
 544     return true;\r
 545 }\r
 546 \r
 547 inline\r
 548 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {\r
 549   \r
 550     if ( SupportData.HasCoreOnly ) return false;\r
 551     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 552     if ( type != "Z" && type != "H" ) return false;\r
 553   \r
 554     // localize the tag data\r
 555     char* pOriginalTagData = (char*)TagData.data();\r
 556     char* pTagData = pOriginalTagData;\r
 557     const unsigned int originalTagDataLength = TagData.size();\r
 558     \r
 559     unsigned int newTagDataLength = 0;\r
 560     unsigned int numBytesParsed = 0;\r
 561     \r
 562     // if tag found, store data in readGroup, return success\r
 563     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {\r
 564         \r
 565         // make sure array is more than big enough\r
 566         char newTagData[originalTagDataLength + value.size()];  \r
 567 \r
 568         // copy original tag data up til desired tag\r
 569         const unsigned int beginningTagDataLength = numBytesParsed;\r
 570         newTagDataLength += beginningTagDataLength;\r
 571         memcpy(newTagData, pOriginalTagData, numBytesParsed);\r
 572       \r
 573         // copy new VALUE in place of current tag data\r
 574         const unsigned int dataLength = strlen(value.c_str());\r
 575         memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );\r
 576         \r
 577         // skip to next tag (if tag for removal is last, return true) \r
 578         const char* pTagStorageType = pTagData - 1;\r
 579         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;\r
 580          \r
 581         // copy everything from current tag (the next one after tag for removal) to end\r
 582         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);\r
 583         const unsigned int endTagOffset      = beginningTagDataLength + dataLength + 1;\r
 584         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;\r
 585         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);\r
 586         \r
 587         // ensure null-terminator\r
 588         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;\r
 589         \r
 590         // save new tag data\r
 591         TagData.assign(newTagData, endTagOffset + endTagDataLength);\r
 592         return true;\r
 593     }\r
 594     \r
 595     // tag not found, attempt AddTag\r
 596     else return AddTag(tag, type, value);\r
 597 }\r
 598 \r
 599 inline\r
 600 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {\r
 601   \r
 602     if ( SupportData.HasCoreOnly ) return false;\r
 603     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 604     if ( type == "f" || type == "Z" || type == "H" ) return false;\r
 605     \r
 606      // localize the tag data\r
 607     char* pOriginalTagData = (char*)TagData.data();\r
 608     char* pTagData = pOriginalTagData;\r
 609     const unsigned int originalTagDataLength = TagData.size();\r
 610     \r
 611     unsigned int newTagDataLength = 0;\r
 612     unsigned int numBytesParsed = 0;\r
 613     \r
 614     // if tag found, store data in readGroup, return success\r
 615     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {\r
 616         \r
 617         // make sure array is more than big enough\r
 618         char newTagData[originalTagDataLength + sizeof(value)];  \r
 619 \r
 620         // copy original tag data up til desired tag\r
 621         const unsigned int beginningTagDataLength = numBytesParsed;\r
 622         newTagDataLength += beginningTagDataLength;\r
 623         memcpy(newTagData, pOriginalTagData, numBytesParsed);\r
 624       \r
 625         // copy new VALUE in place of current tag data\r
 626         union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;\r
 627         un.value = value;\r
 628         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));\r
 629         \r
 630         // skip to next tag (if tag for removal is last, return true) \r
 631         const char* pTagStorageType = pTagData - 1;\r
 632         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;\r
 633          \r
 634         // copy everything from current tag (the next one after tag for removal) to end\r
 635         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);\r
 636         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(unsigned int);\r
 637         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;\r
 638         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);\r
 639         \r
 640         // ensure null-terminator\r
 641         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;\r
 642         \r
 643         // save new tag data\r
 644         TagData.assign(newTagData, endTagOffset + endTagDataLength);\r
 645         return true;\r
 646     }\r
 647     \r
 648     // tag not found, attempt AddTag\r
 649     else return AddTag(tag, type, value);\r
 650 }\r
 651 \r
 652 inline\r
 653 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {\r
 654     return EditTag(tag, type, (const uint32_t&)value);\r
 655 }\r
 656 \r
 657 inline\r
 658 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {\r
 659   \r
 660     if ( SupportData.HasCoreOnly ) return false;\r
 661     if ( tag.size() != 2 || type.size() != 1 ) return false;\r
 662     if ( type == "Z" || type == "H" ) return false;\r
 663     \r
 664      // localize the tag data\r
 665     char* pOriginalTagData = (char*)TagData.data();\r
 666     char* pTagData = pOriginalTagData;\r
 667     const unsigned int originalTagDataLength = TagData.size();\r
 668     \r
 669     unsigned int newTagDataLength = 0;\r
 670     unsigned int numBytesParsed = 0;\r
 671     \r
 672     // if tag found, store data in readGroup, return success\r
 673     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {\r
 674         \r
 675         // make sure array is more than big enough\r
 676         char newTagData[originalTagDataLength + sizeof(value)];  \r
 677 \r
 678         // copy original tag data up til desired tag\r
 679         const unsigned int beginningTagDataLength = numBytesParsed;\r
 680         newTagDataLength += beginningTagDataLength;\r
 681         memcpy(newTagData, pOriginalTagData, numBytesParsed);\r
 682       \r
 683         // copy new VALUE in place of current tag data\r
 684         union { float value; char valueBuffer[sizeof(float)]; } un;\r
 685         un.value = value;\r
 686         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));\r
 687         \r
 688         // skip to next tag (if tag for removal is last, return true) \r
 689         const char* pTagStorageType = pTagData - 1;\r
 690         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;\r
 691          \r
 692         // copy everything from current tag (the next one after tag for removal) to end\r
 693         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);\r
 694         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(float);\r
 695         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;\r
 696         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);\r
 697         \r
 698         // ensure null-terminator\r
 699         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;\r
 700         \r
 701         // save new tag data\r
 702         TagData.assign(newTagData, endTagOffset + endTagDataLength);\r
 703         return true;\r
 704     }\r
 705     \r
 706     // tag not found, attempt AddTag\r
 707     else return AddTag(tag, type, value);\r
 708 }\r
 709 \r
 710 // get "NM" tag data - originally contributed by Aaron Quinlan\r
 711 // stores data in 'editDistance', returns success/fail\r
 712 inline \r
 713 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const { \r
 714     return GetTag("NM", (uint32_t&)editDistance);\r
 715 }\r
 716 \r
 717 // get "RG" tag data\r
 718 // stores data in 'readGroup', returns success/fail\r
 719 inline \r
 720 bool BamAlignment::GetReadGroup(std::string& readGroup) const {\r
 721     return GetTag("RG", readGroup);\r
 722 }\r
 723 \r
 724 inline\r
 725 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {\r
 726 \r
 727     // make sure tag data exists\r
 728     if ( SupportData.HasCoreOnly || TagData.empty() ) \r
 729         return false;\r
 730 \r
 731     // localize the tag data\r
 732     char* pTagData = (char*)TagData.data();\r
 733     const unsigned int tagDataLength = TagData.size();\r
 734     unsigned int numBytesParsed = 0;\r
 735     \r
 736     // if tag found, store data in readGroup, return success\r
 737     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {\r
 738         const unsigned int dataLength = strlen(pTagData);\r
 739         destination.clear();\r
 740         destination.resize(dataLength);\r
 741         memcpy( (char*)destination.data(), pTagData, dataLength );\r
 742         return true;\r
 743     }\r
 744     \r
 745     // tag not found, return failure\r
 746     return false;\r
 747 }\r
 748 \r
 749 inline\r
 750 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {\r
 751   \r
 752     // make sure tag data exists\r
 753     if ( SupportData.HasCoreOnly || TagData.empty() ) \r
 754         return false;\r
 755 \r
 756     // localize the tag data\r
 757     char* pTagData = (char*)TagData.data();\r
 758     const unsigned int tagDataLength = TagData.size();\r
 759     unsigned int numBytesParsed = 0;\r
 760     \r
 761     // if tag found, determine data byte-length, store data in readGroup, return success\r
 762     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {\r
 763         \r
 764         // determine data byte-length\r
 765         const char type = *(pTagData - 1);\r
 766         int destinationLength = 0;\r
 767         switch (type) {\r
 768             // 1 byte data\r
 769             case 'A':\r
 770             case 'c':\r
 771             case 'C':\r
 772                 destinationLength = 1;\r
 773                 break;\r
 774 \r
 775             // 2 byte data\r
 776             case 's':\r
 777             case 'S':\r
 778                 destinationLength = 2;\r
 779                 break;\r
 780 \r
 781             // 4 byte data\r
 782             case 'i':\r
 783             case 'I':\r
 784                 destinationLength = 4;\r
 785                 break;\r
 786 \r
 787             // unsupported type for integer destination (float or var-length strings)\r
 788             case 'f':\r
 789             case 'Z':\r
 790             case 'H':\r
 791                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);\r
 792                 return false;\r
 793 \r
 794             // unknown tag type\r
 795             default:\r
 796                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);\r
 797                 return false;\r
 798         }\r
 799           \r
 800         // store in destination\r
 801         destination = 0;\r
 802         memcpy(&destination, pTagData, destinationLength);\r
 803         return true;\r
 804     }\r
 805     \r
 806     // tag not found, return failure\r
 807     return false;\r
 808 }\r
 809 \r
 810 inline\r
 811 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {\r
 812     return GetTag(tag, (uint32_t&)destination);\r
 813 }\r
 814 \r
 815 inline\r
 816 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {\r
 817   \r
 818     // make sure tag data exists\r
 819     if ( SupportData.HasCoreOnly || TagData.empty() ) \r
 820         return false;\r
 821 \r
 822     // localize the tag data\r
 823     char* pTagData = (char*)TagData.data();\r
 824     const unsigned int tagDataLength = TagData.size();\r
 825     unsigned int numBytesParsed = 0;\r
 826     \r
 827     // if tag found, determine data byte-length, store data in readGroup, return success\r
 828     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {\r
 829         //pTagData += numBytesParsed;\r
 830         \r
 831         // determine data byte-length\r
 832         const char type = *(pTagData - 1);\r
 833         int destinationLength = 0;\r
 834         switch(type) {\r
 835 \r
 836             // 1 byte data\r
 837             case 'A':\r
 838             case 'c':\r
 839             case 'C':\r
 840                 destinationLength = 1;\r
 841                 break;\r
 842 \r
 843             // 2 byte data\r
 844             case 's':\r
 845             case 'S':\r
 846                 destinationLength = 2;\r
 847                 break;\r
 848 \r
 849             // 4 byte data\r
 850             case 'f':\r
 851             case 'i':\r
 852             case 'I':\r
 853                 destinationLength = 4;\r
 854                 break;\r
 855             \r
 856             // unsupported type (var-length strings)\r
 857             case 'Z':\r
 858             case 'H':\r
 859                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);\r
 860                 return false;\r
 861 \r
 862             // unknown tag type\r
 863             default:\r
 864                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);\r
 865                 return false;\r
 866         }\r
 867           \r
 868         // store in destination\r
 869         destination = 0.0;\r
 870         memcpy(&destination, pTagData, destinationLength);\r
 871         return true;\r
 872     }\r
 873     \r
 874     // tag not found, return failure\r
 875     return false;\r
 876 }\r
 877 \r
 878 inline\r
 879 bool BamAlignment::RemoveTag(const std::string& tag) {\r
 880   \r
 881     // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed\r
 882     // also, return false if no data present to remove\r
 883     if ( SupportData.HasCoreOnly || TagData.empty() ) return false;\r
 884   \r
 885     // localize the tag data\r
 886     char* pOriginalTagData = (char*)TagData.data();\r
 887     char* pTagData = pOriginalTagData;\r
 888     const unsigned int originalTagDataLength = TagData.size();\r
 889     unsigned int newTagDataLength = 0;\r
 890     unsigned int numBytesParsed = 0;\r
 891     \r
 892     // if tag found, store data in readGroup, return success\r
 893     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {\r
 894         \r
 895         char newTagData[originalTagDataLength];\r
 896 \r
 897         // copy original tag data up til desired tag\r
 898         pTagData -= 3;\r
 899         numBytesParsed -= 3;\r
 900         const unsigned int beginningTagDataLength = numBytesParsed;\r
 901         newTagDataLength += beginningTagDataLength;\r
 902         memcpy(newTagData, pOriginalTagData, numBytesParsed);\r
 903         \r
 904         // skip to next tag (if tag for removal is last, return true) \r
 905         const char* pTagStorageType = pTagData + 2;\r
 906         pTagData       += 3;\r
 907         numBytesParsed += 3;\r
 908         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;\r
 909          \r
 910         // copy everything from current tag (the next one after tag for removal) to end\r
 911         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);\r
 912         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;\r
 913         memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );\r
 914         \r
 915         // save new tag data\r
 916         TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);\r
 917         return true;\r
 918     }\r
 919     \r
 920     // tag not found, no removal - return failure\r
 921     return false;\r
 922 }\r
 923 \r
 924 inline\r
 925 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {\r
 926 \r
 927     while ( numBytesParsed < tagDataLength ) {\r
 928 \r
 929         const char* pTagType        = pTagData;\r
 930         const char* pTagStorageType = pTagData + 2;\r
 931         pTagData       += 3;\r
 932         numBytesParsed += 3;\r
 933 \r
 934         // check the current tag, return true on match\r
 935         if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 ) \r
 936             return true;\r
 937 \r
 938         // get the storage class and find the next tag\r
 939         if ( *pTagStorageType == '\0' ) return false; \r
 940         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;\r
 941         if ( *pTagData == '\0' ) return false;\r
 942     }\r
 943   \r
 944     // checked all tags, none match\r
 945     return false;\r
 946 }\r
 947 \r
 948 inline\r
 949 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {\r
 950     \r
 951     switch(storageType) {\r
 952 \r
 953         case 'A':\r
 954         case 'c':\r
 955         case 'C':\r
 956             ++numBytesParsed;\r
 957             ++pTagData;\r
 958             break;\r
 959 \r
 960         case 's':\r
 961         case 'S':\r
 962             numBytesParsed += 2;\r
 963             pTagData       += 2;\r
 964             break;\r
 965 \r
 966         case 'f':\r
 967         case 'i':\r
 968         case 'I':\r
 969             numBytesParsed += 4;\r
 970             pTagData       += 4;\r
 971             break;\r
 972 \r
 973         case 'Z':\r
 974         case 'H':\r
 975             while(*pTagData) {\r
 976                 ++numBytesParsed;\r
 977                 ++pTagData;\r
 978             }\r
 979             // increment for null-terminator\r
 980             ++numBytesParsed;\r
 981             ++pTagData;\r
 982             break;\r
 983 \r
 984         default: \r
 985             // error case\r
 986             fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", storageType);\r
 987             return false;\r
 988     }\r
 989     \r
 990     // return success\r
 991     return true;\r
 992 }\r
 993 \r
 994 } // namespace BamTools\r
 995 \r
 996 #endif // BAMAUX_H\r