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: 22 December 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 const char* DNA_LOOKUP = "=ACMGRSVTWYHKDBN";
  24
  25 // default ctor
  26 BamAlignment::BamAlignment(void)
  27     : RefID(-1)
  28     , Position(-1)
  29     , MateRefID(-1)
  30     , MatePosition(-1)
  31     , InsertSize(0)
  32 { }
  33
  34 // copy ctor
  35 BamAlignment::BamAlignment(const BamAlignment& other)
  36     : Name(other.Name)
  37     , Length(other.Length)
  38     , QueryBases(other.QueryBases)
  39     , AlignedBases(other.AlignedBases)
  40     , Qualities(other.Qualities)
  41     , TagData(other.TagData)
  42     , RefID(other.RefID)
  43     , Position(other.Position)
  44     , Bin(other.Bin)
  45     , MapQuality(other.MapQuality)
  46     , AlignmentFlag(other.AlignmentFlag)
  47     , CigarData(other.CigarData)
  48     , MateRefID(other.MateRefID)
  49     , MatePosition(other.MatePosition)
  50     , InsertSize(other.InsertSize)
  51     , SupportData(other.SupportData)
  52 { }
  53
  54 // dtor
  55 BamAlignment::~BamAlignment(void) { }
  56
  57 // Queries against alignment flags
  58 bool BamAlignment::IsDuplicate(void) const         { return ( (AlignmentFlag & DUPLICATE)     != 0 ); }
  59 bool BamAlignment::IsFailedQC(void) const          { return ( (AlignmentFlag & QC_FAILED)     != 0 ); }
  60 bool BamAlignment::IsFirstMate(void) const         { return ( (AlignmentFlag & READ_1)        != 0 ); }
  61 bool BamAlignment::IsMapped(void) const            { return ( (AlignmentFlag & UNMAPPED)      == 0 ); }
  62 bool BamAlignment::IsMateMapped(void) const        { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }
  63 bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }
  64 bool BamAlignment::IsPaired(void) const            { return ( (AlignmentFlag & PAIRED)        != 0 ); }
  65 bool BamAlignment::IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY)     == 0 ); }
  66 bool BamAlignment::IsProperPair(void) const        { return ( (AlignmentFlag & PROPER_PAIR)   != 0 ); }
  67 bool BamAlignment::IsReverseStrand(void) const     { return ( (AlignmentFlag & REVERSE)       != 0 ); }
  68 bool BamAlignment::IsSecondMate(void) const        { return ( (AlignmentFlag & READ_2)        != 0 ); }
  69
  70 // Manipulate alignment flags
  71 void BamAlignment::SetIsDuplicate(bool ok)          { if (ok) AlignmentFlag |= DUPLICATE;     else AlignmentFlag &= ~DUPLICATE; }
  72 void BamAlignment::SetIsFailedQC(bool ok)           { if (ok) AlignmentFlag |= QC_FAILED;     else AlignmentFlag &= ~QC_FAILED; }
  73 void BamAlignment::SetIsFirstMate(bool ok)          { if (ok) AlignmentFlag |= READ_1;        else AlignmentFlag &= ~READ_1; }
  74 void BamAlignment::SetIsMapped(bool ok)             { SetIsUnmapped(!ok); }
  75 void BamAlignment::SetIsMateMapped(bool ok)         { SetIsMateUnmapped(!ok); }
  76 void BamAlignment::SetIsMateUnmapped(bool ok)       { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }
  77 void BamAlignment::SetIsMateReverseStrand(bool ok)  { if (ok) AlignmentFlag |= MATE_REVERSE;  else AlignmentFlag &= ~MATE_REVERSE; }
  78 void BamAlignment::SetIsPaired(bool ok)             { if (ok) AlignmentFlag |= PAIRED;        else AlignmentFlag &= ~PAIRED; }
  79 void BamAlignment::SetIsPrimaryAlignment(bool ok)   { SetIsSecondaryAlignment(!ok); }
  80 void BamAlignment::SetIsProperPair(bool ok)         { if (ok) AlignmentFlag |= PROPER_PAIR;   else AlignmentFlag &= ~PROPER_PAIR; }
  81 void BamAlignment::SetIsReverseStrand(bool ok)      { if (ok) AlignmentFlag |= REVERSE;       else AlignmentFlag &= ~REVERSE; }
  82 void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY;     else AlignmentFlag &= ~SECONDARY; }
  83 void BamAlignment::SetIsSecondMate(bool ok)         { if (ok) AlignmentFlag |= READ_2;        else AlignmentFlag &= ~READ_2; }
  84 void BamAlignment::SetIsUnmapped(bool ok)           { if (ok) AlignmentFlag |= UNMAPPED;      else AlignmentFlag &= ~UNMAPPED; }
  85
  86 // fills out character data
  87 bool BamAlignment::BuildCharData(void) {
  88
  89     // skip if char data already parsed
  90     if ( !SupportData.HasCoreOnly ) return true;
  91
  92     // check system endianness
  93     bool IsBigEndian = BamTools::SystemIsBigEndian();
  94
  95     // calculate character lengths/offsets
  96     const unsigned int dataLength     = SupportData.BlockLength - BAM_CORE_SIZE;
  97     const unsigned int seqDataOffset  = SupportData.QueryNameLength + (SupportData.NumCigarOperations * 4);
  98     const unsigned int qualDataOffset = seqDataOffset + (SupportData.QuerySequenceLength+1)/2;
  99     const unsigned int tagDataOffset  = qualDataOffset + SupportData.QuerySequenceLength;
 100     const unsigned int tagDataLength  = dataLength - tagDataOffset;
 101
 102     // check offsets to see what char data exists
 103     const bool hasSeqData  = ( seqDataOffset  < dataLength );
 104     const bool hasQualData = ( qualDataOffset < dataLength );
 105     const bool hasTagData  = ( tagDataOffset  < dataLength );
 106
 107     // set up char buffers
 108     const char* allCharData = SupportData.AllCharData.data();
 109     const char* seqData     = ( hasSeqData  ? (((const char*)allCharData) + seqDataOffset)  : (const char*)0 );
 110     const char* qualData    = ( hasQualData ? (((const char*)allCharData) + qualDataOffset) : (const char*)0 );
 111           char* tagData     = ( hasTagData  ? (((char*)allCharData) + tagDataOffset)        : (char*)0 );
 112
 113     // store alignment name (relies on null char in name as terminator)
 114     Name.assign((const char*)(allCharData));
 115
 116     // save query sequence
 117     QueryBases.clear();
 118     if ( hasSeqData ) {
 119         QueryBases.reserve(SupportData.QuerySequenceLength);
 120         for (unsigned int i = 0; i < SupportData.QuerySequenceLength; ++i) {
 121             char singleBase = DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
 122             QueryBases.append(1, singleBase);
 123         }
 124     }
 125
 126     // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
 127     Qualities.clear();
 128     if ( hasQualData ) {
 129         Qualities.reserve(SupportData.QuerySequenceLength);
 130         for (unsigned int i = 0; i < SupportData.QuerySequenceLength; ++i) {
 131             char singleQuality = (char)(qualData[i]+33);
 132             Qualities.append(1, singleQuality);
 133         }
 134     }
 135
 136     // clear previous AlignedBases
 137     AlignedBases.clear();
 138
 139     // if QueryBases has data, build AlignedBases using CIGAR data
 140     // otherwise, AlignedBases will remain empty (this case IS allowed)
 141     if ( !QueryBases.empty() ) {
 142
 143         // resize AlignedBases
 144         AlignedBases.reserve(SupportData.QuerySequenceLength);
 145
 146         // iterate over CigarOps
 147         int k = 0;
 148         vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 149         vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 150         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 151
 152             const CigarOp& op = (*cigarIter);
 153             switch(op.Type) {
 154
 155                 // for 'M', 'I' - write bases
 156                 case ('M') :
 157                 case ('I') :
 158                     AlignedBases.append(QueryBases.substr(k, op.Length));
 159                     // fall through
 160
 161                 // for 'S' - soft clip, do not write bases
 162                 // but increment placeholder 'k'
 163                 case ('S') :
 164                     k += op.Length;
 165                     break;
 166
 167                 // for 'D' - write gap character
 168                 case ('D') :
 169                     AlignedBases.append(op.Length, '-');
 170                     break;
 171
 172                 // for 'P' - write padding character
 173                 case ('P') :
 174                     AlignedBases.append( op.Length, '*' );
 175                     break;
 176
 177                 // for 'N' - write N's, skip bases in original query sequence
 178                 case ('N') :
 179                     AlignedBases.append( op.Length, 'N' );
 180                     break;
 181
 182                 // for 'H' - hard clip, do nothing to AlignedBases, move to next op
 183                 case ('H') :
 184                     break;
 185
 186                 // shouldn't get here
 187                 default:
 188                     fprintf(stderr, "ERROR: Invalid Cigar op type\n");
 189                     exit(1);
 190             }
 191         }
 192     }
 193
 194     // save tag data
 195     TagData.clear();
 196     if ( hasTagData ) {
 197         if ( IsBigEndian ) {
 198             int i = 0;
 199             while ( (unsigned int)i < tagDataLength ) {
 200
 201                 i += 2;                                 // skip tagType chars (e.g. "RG", "NM", etc.)
 202                 uint8_t type = toupper(tagData[i]);     // lower & upper case letters have same meaning
 203                 ++i;                                    // skip valueType char (e.g. 'A', 'I', 'Z', etc.)
 204
 205                 switch (type) {
 206
 207                     case('A') :
 208                     case('C') :
 209                         ++i;
 210                         break;
 211
 212                     case('S') :
 213                         SwapEndian_16p(&tagData[i]);
 214                         i += sizeof(uint16_t);
 215                         break;
 216
 217                     case('F') :
 218                     case('I') :
 219                         SwapEndian_32p(&tagData[i]);
 220                         i += sizeof(uint32_t);
 221                         break;
 222
 223                     case('D') :
 224                         SwapEndian_64p(&tagData[i]);
 225                         i += sizeof(uint64_t);
 226                         break;
 227
 228                     case('H') :
 229                     case('Z') :
 230                         while (tagData[i]) { ++i; }
 231                         ++i; // increment one more for null terminator
 232                         break;
 233
 234                     // shouldn't get here
 235                     default :
 236                         fprintf(stderr, "ERROR: Invalid tag value type\n");
 237                         exit(1);
 238                 }
 239             }
 240         }
 241
 242         // store tagData in alignment
 243         TagData.resize(tagDataLength);
 244         memcpy((char*)TagData.data(), tagData, tagDataLength);
 245     }
 246
 247     // clear the core-only flag
 248     SupportData.HasCoreOnly = false;
 249
 250     // return success
 251     return true;
 252 }
 253
 254 // calculates alignment end position, based on starting position and CIGAR operations
 255 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {
 256
 257     // initialize alignment end to starting position
 258     int alignEnd = Position;
 259
 260     // iterate over cigar operations
 261     vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 262     vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 263     for ( ; cigarIter != cigarEnd; ++cigarIter) {
 264         const char cigarType = (*cigarIter).Type;
 265         if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' )
 266             alignEnd += (*cigarIter).Length;
 267         else if ( usePadded && cigarType == 'I' )
 268             alignEnd += (*cigarIter).Length;
 269     }
 270
 271     // adjust for zeroBased, if necessary
 272     if (zeroBased)
 273         return alignEnd - 1;
 274     else
 275         return alignEnd;
 276 }
 277
 278 bool BamAlignment::AddTag(const string& tag, const string& type, const string& value) {
 279
 280     if ( SupportData.HasCoreOnly ) return false;
 281     if ( tag.size() != 2 || type.size() != 1 ) return false;
 282     if ( type != "Z" && type != "H" ) return false;
 283
 284     // localize the tag data
 285     char* pTagData = (char*)TagData.data();
 286     const unsigned int tagDataLength = TagData.size();
 287     unsigned int numBytesParsed = 0;
 288
 289     // if tag already exists, return false
 290     // use EditTag explicitly instead
 291     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 292
 293     // otherwise, copy tag data to temp buffer
 294     string newTag = tag + type + value;
 295     const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
 296     char originalTagData[newTagDataLength];
 297     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 298
 299     // append newTag
 300     strcat(originalTagData + tagDataLength, newTag.data());  // removes original null-term, appends newTag + null-term
 301
 302     // store temp buffer back in TagData
 303     const char* newTagData = (const char*)originalTagData;
 304     TagData.assign(newTagData, newTagDataLength);
 305
 306     // return success
 307     return true;
 308 }
 309
 310 bool BamAlignment::AddTag(const string& tag, const string& type, const uint32_t& value) {
 311
 312     if ( SupportData.HasCoreOnly ) return false;
 313     if ( tag.size() != 2 || type.size() != 1 ) return false;
 314     if ( type == "f" || type == "Z" || type == "H" ) return false;
 315
 316     // localize the tag data
 317     char* pTagData = (char*)TagData.data();
 318     const unsigned int tagDataLength = TagData.size();
 319     unsigned int numBytesParsed = 0;
 320
 321     // if tag already exists, return false
 322     // use EditTag explicitly instead
 323     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 324
 325     // otherwise, convert value to string
 326     union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
 327     un.value = value;
 328
 329     // copy original tag data to temp buffer
 330     string newTag = tag + type;
 331     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
 332     char originalTagData[newTagDataLength];
 333     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 334
 335     // append newTag
 336     strcat(originalTagData + tagDataLength, newTag.data());
 337     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
 338
 339     // store temp buffer back in TagData
 340     const char* newTagData = (const char*)originalTagData;
 341     TagData.assign(newTagData, newTagDataLength);
 342
 343     // return success
 344     return true;
 345 }
 346
 347 bool BamAlignment::AddTag(const string& tag, const string& type, const int32_t& value) {
 348     return AddTag(tag, type, (const uint32_t&)value);
 349 }
 350
 351 bool BamAlignment::AddTag(const string& tag, const string& type, const float& value) {
 352
 353     if ( SupportData.HasCoreOnly ) return false;
 354     if ( tag.size() != 2 || type.size() != 1 ) return false;
 355     if ( type == "Z" || type == "H" ) return false;
 356
 357     // localize the tag data
 358     char* pTagData = (char*)TagData.data();
 359     const unsigned int tagDataLength = TagData.size();
 360     unsigned int numBytesParsed = 0;
 361
 362     // if tag already exists, return false
 363     // use EditTag explicitly instead
 364     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
 365
 366     // otherwise, convert value to string
 367     union { float value; char valueBuffer[sizeof(float)]; } un;
 368     un.value = value;
 369
 370     // copy original tag data to temp buffer
 371     string newTag = tag + type;
 372     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
 373     char originalTagData[newTagDataLength];
 374     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 375
 376     // append newTag
 377     strcat(originalTagData + tagDataLength, newTag.data());
 378     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
 379
 380     // store temp buffer back in TagData
 381     const char* newTagData = (const char*)originalTagData;
 382     TagData.assign(newTagData, newTagDataLength);
 383
 384     // return success
 385     return true;
 386 }
 387
 388 bool BamAlignment::EditTag(const string& tag, const string& type, const string& value) {
 389
 390     if ( SupportData.HasCoreOnly ) return false;
 391     if ( tag.size() != 2 || type.size() != 1 ) return false;
 392     if ( type != "Z" && type != "H" ) return false;
 393
 394     // localize the tag data
 395     char* pOriginalTagData = (char*)TagData.data();
 396     char* pTagData = pOriginalTagData;
 397     const unsigned int originalTagDataLength = TagData.size();
 398
 399     unsigned int newTagDataLength = 0;
 400     unsigned int numBytesParsed = 0;
 401
 402     // if tag found, store data in readGroup, return success
 403     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 404
 405         // make sure array is more than big enough
 406         char newTagData[originalTagDataLength + value.size()];
 407
 408         // copy original tag data up til desired tag
 409         const unsigned int beginningTagDataLength = numBytesParsed;
 410         newTagDataLength += beginningTagDataLength;
 411         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 412
 413         // copy new VALUE in place of current tag data
 414         const unsigned int dataLength = strlen(value.c_str());
 415         memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
 416
 417         // skip to next tag (if tag for removal is last, return true)
 418         const char* pTagStorageType = pTagData - 1;
 419         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 420
 421         // copy everything from current tag (the next one after tag for removal) to end
 422         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 423         const unsigned int endTagOffset      = beginningTagDataLength + dataLength + 1;
 424         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 425         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 426
 427         // ensure null-terminator
 428         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 429
 430         // save new tag data
 431         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 432         return true;
 433     }
 434
 435     // tag not found, attempt AddTag
 436     else return AddTag(tag, type, value);
 437 }
 438
 439 bool BamAlignment::EditTag(const string& tag, const string& type, const uint32_t& value) {
 440
 441     if ( SupportData.HasCoreOnly ) return false;
 442     if ( tag.size() != 2 || type.size() != 1 ) return false;
 443     if ( type == "f" || type == "Z" || type == "H" ) return false;
 444
 445      // localize the tag data
 446     char* pOriginalTagData = (char*)TagData.data();
 447     char* pTagData = pOriginalTagData;
 448     const unsigned int originalTagDataLength = TagData.size();
 449
 450     unsigned int newTagDataLength = 0;
 451     unsigned int numBytesParsed = 0;
 452
 453     // if tag found, store data in readGroup, return success
 454     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 455
 456         // make sure array is more than big enough
 457         char newTagData[originalTagDataLength + sizeof(value)];
 458
 459         // copy original tag data up til desired tag
 460         const unsigned int beginningTagDataLength = numBytesParsed;
 461         newTagDataLength += beginningTagDataLength;
 462         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 463
 464         // copy new VALUE in place of current tag data
 465         union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
 466         un.value = value;
 467         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
 468
 469         // skip to next tag (if tag for removal is last, return true)
 470         const char* pTagStorageType = pTagData - 1;
 471         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 472
 473         // copy everything from current tag (the next one after tag for removal) to end
 474         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 475         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(unsigned int);
 476         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 477         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 478
 479         // ensure null-terminator
 480         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 481
 482         // save new tag data
 483         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 484         return true;
 485     }
 486
 487     // tag not found, attempt AddTag
 488     else return AddTag(tag, type, value);
 489 }
 490
 491 bool BamAlignment::EditTag(const string& tag, const string& type, const int32_t& value) {
 492     return EditTag(tag, type, (const uint32_t&)value);
 493 }
 494
 495 bool BamAlignment::EditTag(const string& tag, const string& type, const float& value) {
 496
 497     if ( SupportData.HasCoreOnly ) return false;
 498     if ( tag.size() != 2 || type.size() != 1 ) return false;
 499     if ( type == "Z" || type == "H" ) return false;
 500
 501      // localize the tag data
 502     char* pOriginalTagData = (char*)TagData.data();
 503     char* pTagData = pOriginalTagData;
 504     const unsigned int originalTagDataLength = TagData.size();
 505
 506     unsigned int newTagDataLength = 0;
 507     unsigned int numBytesParsed = 0;
 508
 509     // if tag found, store data in readGroup, return success
 510     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 511
 512         // make sure array is more than big enough
 513         char newTagData[originalTagDataLength + sizeof(value)];
 514
 515         // copy original tag data up til desired tag
 516         const unsigned int beginningTagDataLength = numBytesParsed;
 517         newTagDataLength += beginningTagDataLength;
 518         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 519
 520         // copy new VALUE in place of current tag data
 521         union { float value; char valueBuffer[sizeof(float)]; } un;
 522         un.value = value;
 523         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
 524
 525         // skip to next tag (if tag for removal is last, return true)
 526         const char* pTagStorageType = pTagData - 1;
 527         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 528
 529         // copy everything from current tag (the next one after tag for removal) to end
 530         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 531         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(float);
 532         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 533         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 534
 535         // ensure null-terminator
 536         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 537
 538         // save new tag data
 539         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 540         return true;
 541     }
 542
 543     // tag not found, attempt AddTag
 544     else return AddTag(tag, type, value);
 545 }
 546
 547 // get "NM" tag data - originally contributed by Aaron Quinlan
 548 // stores data in 'editDistance', returns success/fail
 549 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
 550     return GetTag("NM", (uint32_t&)editDistance);
 551 }
 552
 553 // get "RG" tag data
 554 // stores data in 'readGroup', returns success/fail
 555 bool BamAlignment::GetReadGroup(string& readGroup) const {
 556     return GetTag("RG", readGroup);
 557 }
 558
 559 bool BamAlignment::GetTag(const string& tag, string& destination) const {
 560
 561     // make sure tag data exists
 562     if ( SupportData.HasCoreOnly || TagData.empty() )
 563         return false;
 564
 565     // localize the tag data
 566     char* pTagData = (char*)TagData.data();
 567     const unsigned int tagDataLength = TagData.size();
 568     unsigned int numBytesParsed = 0;
 569
 570     // if tag found, store data in readGroup, return success
 571     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 572         const unsigned int dataLength = strlen(pTagData);
 573         destination.clear();
 574         destination.resize(dataLength);
 575         memcpy( (char*)destination.data(), pTagData, dataLength );
 576         return true;
 577     }
 578
 579     // tag not found, return failure
 580     return false;
 581 }
 582
 583 bool BamAlignment::GetTag(const string& tag, uint32_t& destination) const {
 584
 585     // make sure tag data exists
 586     if ( SupportData.HasCoreOnly || TagData.empty() )
 587         return false;
 588
 589     // localize the tag data
 590     char* pTagData = (char*)TagData.data();
 591     const unsigned int tagDataLength = TagData.size();
 592     unsigned int numBytesParsed = 0;
 593
 594     // if tag found, determine data byte-length, store data in readGroup, return success
 595     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 596
 597         // determine data byte-length
 598         const char type = *(pTagData - 1);
 599         int destinationLength = 0;
 600         switch (type) {
 601
 602             // 1 byte data
 603             case 'A':
 604             case 'c':
 605             case 'C':
 606                 destinationLength = 1;
 607                 break;
 608
 609             // 2 byte data
 610             case 's':
 611             case 'S':
 612                 destinationLength = 2;
 613                 break;
 614
 615             // 4 byte data
 616             case 'i':
 617             case 'I':
 618                 destinationLength = 4;
 619                 break;
 620
 621             // unsupported type for integer destination (float or var-length strings)
 622             case 'f':
 623             case 'Z':
 624             case 'H':
 625                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
 626                 return false;
 627
 628             // unknown tag type
 629             default:
 630                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
 631                 return false;
 632         }
 633
 634         // store in destination
 635         destination = 0;
 636         memcpy(&destination, pTagData, destinationLength);
 637         return true;
 638     }
 639
 640     // tag not found, return failure
 641     return false;
 642 }
 643
 644 bool BamAlignment::GetTag(const string& tag, int32_t& destination) const {
 645     return GetTag(tag, (uint32_t&)destination);
 646 }
 647
 648 bool BamAlignment::GetTag(const string& tag, float& destination) const {
 649
 650     // make sure tag data exists
 651     if ( SupportData.HasCoreOnly || TagData.empty() )
 652         return false;
 653
 654     // localize the tag data
 655     char* pTagData = (char*)TagData.data();
 656     const unsigned int tagDataLength = TagData.size();
 657     unsigned int numBytesParsed = 0;
 658
 659     // if tag found, determine data byte-length, store data in readGroup, return success
 660     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 661
 662         // determine data byte-length
 663         const char type = *(pTagData - 1);
 664         int destinationLength = 0;
 665         switch(type) {
 666
 667             // 1 byte data
 668             case 'A':
 669             case 'c':
 670             case 'C':
 671                 destinationLength = 1;
 672                 break;
 673
 674             // 2 byte data
 675             case 's':
 676             case 'S':
 677                 destinationLength = 2;
 678                 break;
 679
 680             // 4 byte data
 681             case 'f':
 682             case 'i':
 683             case 'I':
 684                 destinationLength = 4;
 685                 break;
 686
 687             // unsupported type (var-length strings)
 688             case 'Z':
 689             case 'H':
 690                 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
 691                 return false;
 692
 693             // unknown tag type
 694             default:
 695                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
 696                 return false;
 697         }
 698
 699         // store in destination
 700         destination = 0.0;
 701         memcpy(&destination, pTagData, destinationLength);
 702         return true;
 703     }
 704
 705     // tag not found, return failure
 706     return false;
 707 }
 708
 709 bool BamAlignment::GetTagType(const string& tag, char& type) const {
 710
 711     // make sure tag data exists
 712     if ( SupportData.HasCoreOnly || TagData.empty() )
 713         return false;
 714
 715     // localize the tag data
 716     char* pTagData = (char*)TagData.data();
 717     const unsigned int tagDataLength = TagData.size();
 718     unsigned int numBytesParsed = 0;
 719
 720     // lookup tag
 721     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 722
 723         // retrieve tag type code
 724         type = *(pTagData - 1);
 725
 726         // validate that type is a proper BAM tag type
 727         switch(type) {
 728             case 'A':
 729             case 'c':
 730             case 'C':
 731             case 's':
 732             case 'S':
 733             case 'f':
 734             case 'i':
 735             case 'I':
 736             case 'Z':
 737             case 'H':
 738                 return true;
 739
 740             // unknown tag type
 741             default:
 742                 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
 743                 return false;
 744         }
 745     }
 746
 747     // tag not found, return failure
 748     return false;
 749 }
 750
 751 bool BamAlignment::RemoveTag(const string& tag) {
 752
 753     // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
 754     // also, return false if no data present to remove
 755     if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
 756
 757     // localize the tag data
 758     char* pOriginalTagData = (char*)TagData.data();
 759     char* pTagData = pOriginalTagData;
 760     const unsigned int originalTagDataLength = TagData.size();
 761     unsigned int newTagDataLength = 0;
 762     unsigned int numBytesParsed = 0;
 763
 764     // if tag found, store data in readGroup, return success
 765     if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 766
 767         char newTagData[originalTagDataLength];
 768
 769         // copy original tag data up til desired tag
 770         pTagData -= 3;
 771         numBytesParsed -= 3;
 772         const unsigned int beginningTagDataLength = numBytesParsed;
 773         newTagDataLength += beginningTagDataLength;
 774         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 775
 776         // skip to next tag (if tag for removal is last, return true)
 777         const char* pTagStorageType = pTagData + 2;
 778         pTagData       += 3;
 779         numBytesParsed += 3;
 780         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
 781
 782         // copy everything from current tag (the next one after tag for removal) to end
 783         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 784         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 785         memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
 786
 787         // save new tag data
 788         TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
 789         return true;
 790     }
 791
 792     // tag not found, no removal - return failure
 793     return false;
 794 }
 795
 796 bool BamAlignment::FindTag(const string& tag,
 797                            char* &pTagData,
 798                            const unsigned int& tagDataLength,
 799                            unsigned int& numBytesParsed)
 800 {
 801
 802     while ( numBytesParsed < tagDataLength ) {
 803
 804         const char* pTagType        = pTagData;
 805         const char* pTagStorageType = pTagData + 2;
 806         pTagData       += 3;
 807         numBytesParsed += 3;
 808
 809         // check the current tag, return true on match
 810         if ( strncmp(pTagType, tag.c_str(), 2) == 0 )
 811             return true;
 812
 813         // get the storage class and find the next tag
 814         if ( *pTagStorageType == '\0' ) return false;
 815         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
 816         if ( *pTagData == '\0' ) return false;
 817     }
 818
 819     // checked all tags, none match
 820     return false;
 821 }
 822
 823 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
 824
 825     switch(storageType) {
 826
 827         case 'A':
 828         case 'c':
 829         case 'C':
 830             ++numBytesParsed;
 831             ++pTagData;
 832             break;
 833
 834         case 's':
 835         case 'S':
 836             numBytesParsed += 2;
 837             pTagData       += 2;
 838             break;
 839
 840         case 'f':
 841         case 'i':
 842         case 'I':
 843             numBytesParsed += 4;
 844             pTagData       += 4;
 845             break;
 846
 847         case 'Z':
 848         case 'H':
 849             while(*pTagData) {
 850                 ++numBytesParsed;
 851                 ++pTagData;
 852             }
 853             // increment for null-terminator
 854             ++numBytesParsed;
 855             ++pTagData;
 856             break;
 857
 858         default:
 859             // error case
 860             fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
 861             return false;
 862     }
 863
 864     // return success
 865     return true;
 866 }