src/api/BamAlignment.cpp

   1 // ***************************************************************************
   2 // BamAlignment.cpp (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 18 September 2010 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides the BamAlignment data structure
   9 // ***************************************************************************
  10
  11 #include <cctype>
  12 #include <cstdio>
  13 #include <cstdlib>
  14 #include <cstring>
  15 #include <exception>
  16 #include <map>
  17 #include <utility>
  18 #include "BamAlignment.h"
  19 using namespace BamTools;
  20
  21 // default ctor
  22 BamAlignment::BamAlignment(void)
  23     : RefID(-1)
  24     , Position(-1)
  25     , MateRefID(-1)
  26     , MatePosition(-1)
  27     , InsertSize(0)
  28 { }
  29
  30 // copy ctor
  31 BamAlignment::BamAlignment(const BamAlignment& other)
  32     : Name(other.Name)
  33     , Length(other.Length)
  34     , QueryBases(other.QueryBases)
  35     , AlignedBases(other.AlignedBases)
  36     , Qualities(other.Qualities)
  37     , TagData(other.TagData)
  38     , RefID(other.RefID)
  39     , Position(other.Position)
  40     , Bin(other.Bin)
  41     , MapQuality(other.MapQuality)
  42     , AlignmentFlag(other.AlignmentFlag)
  43     , CigarData(other.CigarData)
  44     , MateRefID(other.MateRefID)
  45     , MatePosition(other.MatePosition)
  46     , InsertSize(other.InsertSize)
  47     , SupportData(other.SupportData)
  48 { }
  49
  50 // dtor
  51 BamAlignment::~BamAlignment(void) { }
  52
  53 // Queries against alignment flags
  54 bool BamAlignment::IsDuplicate(void) const         { return ( (AlignmentFlag & DUPLICATE)     != 0 ); }
  55 bool BamAlignment::IsFailedQC(void) const          { return ( (AlignmentFlag & QC_FAILED)     != 0 ); }
  56 bool BamAlignment::IsFirstMate(void) const         { return ( (AlignmentFlag & READ_1)        != 0 ); }
  57 bool BamAlignment::IsMapped(void) const            { return ( (AlignmentFlag & UNMAPPED)      == 0 ); }
  58 bool BamAlignment::IsMateMapped(void) const        { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }
  59 bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }
  60 bool BamAlignment::IsPaired(void) const            { return ( (AlignmentFlag & PAIRED)        != 0 ); }
  61 bool BamAlignment::IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY)     == 0 ); }
  62 bool BamAlignment::IsProperPair(void) const        { return ( (AlignmentFlag & PROPER_PAIR)   != 0 ); }
  63 bool BamAlignment::IsReverseStrand(void) const     { return ( (AlignmentFlag & REVERSE)       != 0 ); }
  64 bool BamAlignment::IsSecondMate(void) const        { return ( (AlignmentFlag & READ_2)        != 0 ); }
  65
  66 // Manipulate alignment flags
  67 void BamAlignment::SetIsDuplicate(bool ok)          { if (ok) AlignmentFlag |= DUPLICATE;     else AlignmentFlag &= ~DUPLICATE; }
  68 void BamAlignment::SetIsFailedQC(bool ok)           { if (ok) AlignmentFlag |= QC_FAILED;     else AlignmentFlag &= ~QC_FAILED; }
  69 void BamAlignment::SetIsFirstMate(bool ok)          { if (ok) AlignmentFlag |= READ_1;        else AlignmentFlag &= ~READ_1; }
  70 void BamAlignment::SetIsMateUnmapped(bool ok)       { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }
  71 void BamAlignment::SetIsMateReverseStrand(bool ok)  { if (ok) AlignmentFlag |= MATE_REVERSE;  else AlignmentFlag &= ~MATE_REVERSE; }
  72 void BamAlignment::SetIsPaired(bool ok)             { if (ok) AlignmentFlag |= PAIRED;        else AlignmentFlag &= ~PAIRED; }
  73 void BamAlignment::SetIsProperPair(bool ok)         { if (ok) AlignmentFlag |= PROPER_PAIR;   else AlignmentFlag &= ~PROPER_PAIR; }
  74 void BamAlignment::SetIsReverseStrand(bool ok)      { if (ok) AlignmentFlag |= REVERSE;       else AlignmentFlag &= ~REVERSE; }
  75 void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY;     else AlignmentFlag &= ~SECONDARY; }
  76 void BamAlignment::SetIsSecondMate(bool ok)         { if (ok) AlignmentFlag |= READ_2;        else AlignmentFlag &= ~READ_2; }
  77 void BamAlignment::SetIsUnmapped(bool ok)           { if (ok) AlignmentFlag |= UNMAPPED;      else AlignmentFlag &= ~UNMAPPED; }
  78
  79 // calculates alignment end position, based on starting position and CIGAR operations
  80 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {
  81
  82     // initialize alignment end to starting position
  83     int alignEnd = Position;
  84
  85     // iterate over cigar operations
  86     std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
  87     std::vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
  88     for ( ; cigarIter != cigarEnd; ++cigarIter) {
  89         const char cigarType = (*cigarIter).Type;
  90         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' )
  91             alignEnd += (*cigarIter).Length;
  92         else if ( usePadded && cigarType == 'I' )
  93             alignEnd += (*cigarIter).Length;
  94     }
  95
  96     // adjust for zeroBased, if necessary
  97     if (zeroBased)
  98         return alignEnd - 1;
  99     else
 100         return alignEnd;
 101 }
 102
 103 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
 104
 105     if ( SupportData.HasCoreOnly ) return false;
 106     if ( tag.size() != 2 || type.size() != 1 ) return false;
 107     if ( type != "Z" && type != "H" ) return false;
 108
 109     // localize the tag data
 110     char* pTagData = (char*)TagData.data();
 111     const unsigned int tagDataLength = TagData.size();
 112     unsigned int numBytesParsed = 0;
 113
 114     // if tag already exists, return false
 115     // use EditTag explicitly instead
 116     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 117
 118     // otherwise, copy tag data to temp buffer
 119     std::string newTag = tag + type + value;
 120     const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
 121     char originalTagData[newTagDataLength];
 122     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 123
 124     // append newTag
 125     strcat(originalTagData + tagDataLength, newTag.data());  // removes original null-term, appends newTag + null-term
 126
 127     // store temp buffer back in TagData
 128     const char* newTagData = (const char*)originalTagData;
 129     TagData.assign(newTagData, newTagDataLength);
 130
 131     // return success
 132     return true;
 133 }
 134
 135 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
 136
 137     if ( SupportData.HasCoreOnly ) return false;
 138     if ( tag.size() != 2 || type.size() != 1 ) return false;
 139     if ( type == "f" || type == "Z" || type == "H" ) return false;
 140
 141     // localize the tag data
 142     char* pTagData = (char*)TagData.data();
 143     const unsigned int tagDataLength = TagData.size();
 144     unsigned int numBytesParsed = 0;
 145
 146     // if tag already exists, return false
 147     // use EditTag explicitly instead
 148     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 149
 150     // otherwise, convert value to string
 151     union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
 152     un.value = value;
 153
 154     // copy original tag data to temp buffer
 155     std::string newTag = tag + type;
 156     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
 157     char originalTagData[newTagDataLength];
 158     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 159
 160     // append newTag
 161     strcat(originalTagData + tagDataLength, newTag.data());
 162     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
 163
 164     // store temp buffer back in TagData
 165     const char* newTagData = (const char*)originalTagData;
 166     TagData.assign(newTagData, newTagDataLength);
 167
 168     // return success
 169     return true;
 170 }
 171
 172 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
 173     return AddTag(tag, type, (const uint32_t&)value);
 174 }
 175
 176 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
 177
 178     if ( SupportData.HasCoreOnly ) return false;
 179     if ( tag.size() != 2 || type.size() != 1 ) return false;
 180     if ( type == "Z" || type == "H" ) return false;
 181
 182     // localize the tag data
 183     char* pTagData = (char*)TagData.data();
 184     const unsigned int tagDataLength = TagData.size();
 185     unsigned int numBytesParsed = 0;
 186
 187     // if tag already exists, return false
 188     // use EditTag explicitly instead
 189     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 190
 191     // otherwise, convert value to string
 192     union { float value; char valueBuffer[sizeof(float)]; } un;
 193     un.value = value;
 194
 195     // copy original tag data to temp buffer
 196     std::string newTag = tag + type;
 197     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
 198     char originalTagData[newTagDataLength];
 199     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 200
 201     // append newTag
 202     strcat(originalTagData + tagDataLength, newTag.data());
 203     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
 204
 205     // store temp buffer back in TagData
 206     const char* newTagData = (const char*)originalTagData;
 207     TagData.assign(newTagData, newTagDataLength);
 208
 209     // return success
 210     return true;
 211 }
 212
 213 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
 214
 215     if ( SupportData.HasCoreOnly ) return false;
 216     if ( tag.size() != 2 || type.size() != 1 ) return false;
 217     if ( type != "Z" && type != "H" ) return false;
 218
 219     // localize the tag data
 220     char* pOriginalTagData = (char*)TagData.data();
 221     char* pTagData = pOriginalTagData;
 222     const unsigned int originalTagDataLength = TagData.size();
 223
 224     unsigned int newTagDataLength = 0;
 225     unsigned int numBytesParsed = 0;
 226
 227     // if tag found, store data in readGroup, return success
 228     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 229
 230         // make sure array is more than big enough
 231         char newTagData[originalTagDataLength + value.size()];
 232
 233         // copy original tag data up til desired tag
 234         const unsigned int beginningTagDataLength = numBytesParsed;
 235         newTagDataLength += beginningTagDataLength;
 236         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 237
 238         // copy new VALUE in place of current tag data
 239         const unsigned int dataLength = strlen(value.c_str());
 240         memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
 241
 242         // skip to next tag (if tag for removal is last, return true)
 243         const char* pTagStorageType = pTagData - 1;
 244         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 245
 246         // copy everything from current tag (the next one after tag for removal) to end
 247         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 248         const unsigned int endTagOffset      = beginningTagDataLength + dataLength + 1;
 249         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 250         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 251
 252         // ensure null-terminator
 253         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 254
 255         // save new tag data
 256         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 257         return true;
 258     }
 259
 260     // tag not found, attempt AddTag
 261     else return AddTag(tag, type, value);
 262 }
 263
 264 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
 265
 266     if ( SupportData.HasCoreOnly ) return false;
 267     if ( tag.size() != 2 || type.size() != 1 ) return false;
 268     if ( type == "f" || type == "Z" || type == "H" ) return false;
 269
 270      // localize the tag data
 271     char* pOriginalTagData = (char*)TagData.data();
 272     char* pTagData = pOriginalTagData;
 273     const unsigned int originalTagDataLength = TagData.size();
 274
 275     unsigned int newTagDataLength = 0;
 276     unsigned int numBytesParsed = 0;
 277
 278     // if tag found, store data in readGroup, return success
 279     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 280
 281         // make sure array is more than big enough
 282         char newTagData[originalTagDataLength + sizeof(value)];
 283
 284         // copy original tag data up til desired tag
 285         const unsigned int beginningTagDataLength = numBytesParsed;
 286         newTagDataLength += beginningTagDataLength;
 287         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 288
 289         // copy new VALUE in place of current tag data
 290         union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
 291         un.value = value;
 292         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
 293
 294         // skip to next tag (if tag for removal is last, return true)
 295         const char* pTagStorageType = pTagData - 1;
 296         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 297
 298         // copy everything from current tag (the next one after tag for removal) to end
 299         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 300         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(unsigned int);
 301         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 302         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 303
 304         // ensure null-terminator
 305         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 306
 307         // save new tag data
 308         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 309         return true;
 310     }
 311
 312     // tag not found, attempt AddTag
 313     else return AddTag(tag, type, value);
 314 }
 315
 316 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
 317     return EditTag(tag, type, (const uint32_t&)value);
 318 }
 319
 320 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
 321
 322     if ( SupportData.HasCoreOnly ) return false;
 323     if ( tag.size() != 2 || type.size() != 1 ) return false;
 324     if ( type == "Z" || type == "H" ) return false;
 325
 326      // localize the tag data
 327     char* pOriginalTagData = (char*)TagData.data();
 328     char* pTagData = pOriginalTagData;
 329     const unsigned int originalTagDataLength = TagData.size();
 330
 331     unsigned int newTagDataLength = 0;
 332     unsigned int numBytesParsed = 0;
 333
 334     // if tag found, store data in readGroup, return success
 335     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 336
 337         // make sure array is more than big enough
 338         char newTagData[originalTagDataLength + sizeof(value)];
 339
 340         // copy original tag data up til desired tag
 341         const unsigned int beginningTagDataLength = numBytesParsed;
 342         newTagDataLength += beginningTagDataLength;
 343         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 344
 345         // copy new VALUE in place of current tag data
 346         union { float value; char valueBuffer[sizeof(float)]; } un;
 347         un.value = value;
 348         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
 349
 350         // skip to next tag (if tag for removal is last, return true)
 351         const char* pTagStorageType = pTagData - 1;
 352         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 353
 354         // copy everything from current tag (the next one after tag for removal) to end
 355         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 356         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(float);
 357         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 358         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 359
 360         // ensure null-terminator
 361         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 362
 363         // save new tag data
 364         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 365         return true;
 366     }
 367
 368     // tag not found, attempt AddTag
 369     else return AddTag(tag, type, value);
 370 }
 371
 372 // get "NM" tag data - originally contributed by Aaron Quinlan
 373 // stores data in 'editDistance', returns success/fail
 374 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
 375     return GetTag("NM", (uint32_t&)editDistance);
 376 }
 377
 378 // get "RG" tag data
 379 // stores data in 'readGroup', returns success/fail
 380 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
 381     return GetTag("RG", readGroup);
 382 }
 383
 384 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
 385
 386     // make sure tag data exists
 387     if ( SupportData.HasCoreOnly || TagData.empty() )
 388         return false;
 389
 390     // localize the tag data
 391     char* pTagData = (char*)TagData.data();
 392     const unsigned int tagDataLength = TagData.size();
 393     unsigned int numBytesParsed = 0;
 394
 395     // if tag found, store data in readGroup, return success
 396     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 397         const unsigned int dataLength = strlen(pTagData);
 398         destination.clear();
 399         destination.resize(dataLength);
 400         memcpy( (char*)destination.data(), pTagData, dataLength );
 401         return true;
 402     }
 403
 404     // tag not found, return failure
 405     return false;
 406 }
 407
 408 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
 409
 410     // make sure tag data exists
 411     if ( SupportData.HasCoreOnly || TagData.empty() )
 412         return false;
 413
 414     // localize the tag data
 415     char* pTagData = (char*)TagData.data();
 416     const unsigned int tagDataLength = TagData.size();
 417     unsigned int numBytesParsed = 0;
 418
 419     // if tag found, determine data byte-length, store data in readGroup, return success
 420     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 421
 422         // determine data byte-length
 423         const char type = *(pTagData - 1);
 424         int destinationLength = 0;
 425         switch (type) {
 426             // 1 byte data
 427             case 'A':
 428             case 'c':
 429             case 'C':
 430                 destinationLength = 1;
 431                 break;
 432
 433             // 2 byte data
 434             case 's':
 435             case 'S':
 436                 destinationLength = 2;
 437                 break;
 438
 439             // 4 byte data
 440             case 'i':
 441             case 'I':
 442                 destinationLength = 4;
 443                 break;
 444
 445             // unsupported type for integer destination (float or var-length strings)
 446             case 'f':
 447             case 'Z':
 448             case 'H':
 449                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
 450                 return false;
 451
 452             // unknown tag type
 453             default:
 454                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
 455                 return false;
 456         }
 457
 458         // store in destination
 459         destination = 0;
 460         memcpy(&destination, pTagData, destinationLength);
 461         return true;
 462     }
 463
 464     // tag not found, return failure
 465     return false;
 466 }
 467
 468 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
 469     return GetTag(tag, (uint32_t&)destination);
 470 }
 471
 472 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
 473
 474     // make sure tag data exists
 475     if ( SupportData.HasCoreOnly || TagData.empty() )
 476         return false;
 477
 478     // localize the tag data
 479     char* pTagData = (char*)TagData.data();
 480     const unsigned int tagDataLength = TagData.size();
 481     unsigned int numBytesParsed = 0;
 482
 483     // if tag found, determine data byte-length, store data in readGroup, return success
 484     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 485         //pTagData += numBytesParsed;
 486
 487         // determine data byte-length
 488         const char type = *(pTagData - 1);
 489         int destinationLength = 0;
 490         switch(type) {
 491
 492             // 1 byte data
 493             case 'A':
 494             case 'c':
 495             case 'C':
 496                 destinationLength = 1;
 497                 break;
 498
 499             // 2 byte data
 500             case 's':
 501             case 'S':
 502                 destinationLength = 2;
 503                 break;
 504
 505             // 4 byte data
 506             case 'f':
 507             case 'i':
 508             case 'I':
 509                 destinationLength = 4;
 510                 break;
 511
 512             // unsupported type (var-length strings)
 513             case 'Z':
 514             case 'H':
 515                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
 516                 return false;
 517
 518             // unknown tag type
 519             default:
 520                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
 521                 return false;
 522         }
 523
 524         // store in destination
 525         destination = 0.0;
 526         memcpy(&destination, pTagData, destinationLength);
 527         return true;
 528     }
 529
 530     // tag not found, return failure
 531     return false;
 532 }
 533
 534 bool BamAlignment::RemoveTag(const std::string& tag) {
 535
 536     // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
 537     // also, return false if no data present to remove
 538     if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
 539
 540     // localize the tag data
 541     char* pOriginalTagData = (char*)TagData.data();
 542     char* pTagData = pOriginalTagData;
 543     const unsigned int originalTagDataLength = TagData.size();
 544     unsigned int newTagDataLength = 0;
 545     unsigned int numBytesParsed = 0;
 546
 547     // if tag found, store data in readGroup, return success
 548     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 549
 550         char newTagData[originalTagDataLength];
 551
 552         // copy original tag data up til desired tag
 553         pTagData -= 3;
 554         numBytesParsed -= 3;
 555         const unsigned int beginningTagDataLength = numBytesParsed;
 556         newTagDataLength += beginningTagDataLength;
 557         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 558
 559         // skip to next tag (if tag for removal is last, return true)
 560         const char* pTagStorageType = pTagData + 2;
 561         pTagData       += 3;
 562         numBytesParsed += 3;
 563         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 564
 565         // copy everything from current tag (the next one after tag for removal) to end
 566         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 567         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 568         memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
 569
 570         // save new tag data
 571         TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
 572         return true;
 573     }
 574
 575     // tag not found, no removal - return failure
 576     return false;
 577 }
 578
 579 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {
 580
 581     while ( numBytesParsed < tagDataLength ) {
 582
 583         const char* pTagType        = pTagData;
 584         const char* pTagStorageType = pTagData + 2;
 585         pTagData       += 3;
 586         numBytesParsed += 3;
 587
 588         // check the current tag, return true on match
 589         if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 )
 590             return true;
 591
 592         // get the storage class and find the next tag
 593         if ( *pTagStorageType == '\0' ) return false;
 594         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
 595         if ( *pTagData == '\0' ) return false;
 596     }
 597
 598     // checked all tags, none match
 599     return false;
 600 }
 601
 602 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
 603
 604     switch(storageType) {
 605
 606         case 'A':
 607         case 'c':
 608         case 'C':
 609             ++numBytesParsed;
 610             ++pTagData;
 611             break;
 612
 613         case 's':
 614         case 'S':
 615             numBytesParsed += 2;
 616             pTagData       += 2;
 617             break;
 618
 619         case 'f':
 620         case 'i':
 621         case 'I':
 622             numBytesParsed += 4;
 623             pTagData       += 4;
 624             break;
 625
 626         case 'Z':
 627         case 'H':
 628             while(*pTagData) {
 629                 ++numBytesParsed;
 630                 ++pTagData;
 631             }
 632             // increment for null-terminator
 633             ++numBytesParsed;
 634             ++pTagData;
 635             break;
 636
 637         default:
 638             // error case
 639             fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
 640             return false;
 641     }
 642
 643     // return success
 644     return true;
 645 }