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