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