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: 21 March 2011 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides the BamAlignment data structure
   9 // ***************************************************************************
  10
  11 #include <api/BamAlignment.h>
  12 #include <api/BamConstants.h>
  13 using namespace BamTools;
  14
  15 #include <cctype>
  16 #include <cstdio>
  17 #include <cstdlib>
  18 #include <cstring>
  19 #include <exception>
  20 #include <map>
  21 #include <utility>
  22 using namespace std;
  23
  24 // internal utility methods
  25 namespace BamTools {
  26 namespace Internal {
  27
  28 /*! \fn bool IsValidSize(const string& tag, const string& type)
  29     \internal
  30
  31     Checks that tag name & type strings are expected sizes.
  32     \a tag  should have length
  33     \a type should have length 1
  34
  35     \param tag  BAM tag name
  36     \param type BAM tag type-code
  37
  38     \return \c true if both \a tag and \a type are correct sizes
  39 */
  40 bool IsValidSize(const string& tag, const string& type) {
  41     return (tag.size()  == Constants::BAM_TAG_TAGSIZE) &&
  42            (type.size() == Constants::BAM_TAG_TYPESIZE);
  43 }
  44
  45 /*! \fn bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed)
  46     \internal
  47
  48     Moves to next available tag in tag data string
  49
  50     \param storageType    BAM tag type-code that determines how far to move cursor
  51     \param pTagData       pointer to current position (cursor) in tag string
  52     \param numBytesParsed report of how many bytes were parsed (cumulatively)
  53
  54     \return \c if storageType was a recognized BAM tag type
  55     \post \a pTagData will point to the byte where the next tag data begins.
  56           \a numBytesParsed will correspond to the cursor's position in the full TagData string.
  57 */
  58 bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
  59
  60     switch (storageType) {
  61
  62         case (Constants::BAM_TAG_TYPE_ASCII) :
  63         case (Constants::BAM_TAG_TYPE_INT8)  :
  64         case (Constants::BAM_TAG_TYPE_UINT8) :
  65             ++numBytesParsed;
  66             ++pTagData;
  67             break;
  68
  69         case (Constants::BAM_TAG_TYPE_INT16)  :
  70         case (Constants::BAM_TAG_TYPE_UINT16) :
  71             numBytesParsed += 2;
  72             pTagData       += 2;
  73             break;
  74
  75         case (Constants::BAM_TAG_TYPE_FLOAT)  :
  76         case (Constants::BAM_TAG_TYPE_INT32)  :
  77         case (Constants::BAM_TAG_TYPE_UINT32) :
  78             numBytesParsed += 4;
  79             pTagData       += 4;
  80             break;
  81
  82         case (Constants::BAM_TAG_TYPE_STRING) :
  83         case (Constants::BAM_TAG_TYPE_HEX)    :
  84             while(*pTagData) {
  85                 ++numBytesParsed;
  86                 ++pTagData;
  87             }
  88             // increment for null-terminator
  89             ++numBytesParsed;
  90             ++pTagData;
  91             break;
  92
  93         default:
  94             // error case
  95             fprintf(stderr, "BamAlignment ERROR: unknown tag type encountered: [%c]\n", storageType);
  96             return false;
  97     }
  98
  99     // return success
 100     return true;
 101 }
 102
 103 /*! \fn bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed)
 104     \internal
 105
 106     Searches for requested tag in BAM tag data.
 107
 108     \param tag            requested 2-character tag name
 109     \param pTagData       pointer to current position in BamAlignment::TagData
 110     \param tagDataLength  length of BamAlignment::TagData
 111     \param numBytesParsed number of bytes parsed so far
 112
 113     \return \c true if found
 114
 115     \post If \a tag is found, \a pTagData will point to the byte where the tag data begins.
 116           \a numBytesParsed will correspond to the position in the full TagData string.
 117
 118 */
 119 bool FindTag(const std::string& tag,
 120              char* &pTagData,
 121              const unsigned int& tagDataLength,
 122              unsigned int& numBytesParsed)
 123 {
 124
 125     while ( numBytesParsed < tagDataLength ) {
 126
 127         const char* pTagType        = pTagData;
 128         const char* pTagStorageType = pTagData + 2;
 129         pTagData       += 3;
 130         numBytesParsed += 3;
 131
 132         // check the current tag, return true on match
 133         if ( strncmp(pTagType, tag.c_str(), 2) == 0 )
 134             return true;
 135
 136         // get the storage class and find the next tag
 137         if ( *pTagStorageType == '\0' ) return false;
 138         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
 139         if ( *pTagData == '\0' ) return false;
 140     }
 141
 142     // checked all tags, none match
 143     return false;
 144 }
 145
 146 } // namespace Internal
 147 } // namespace BamTools
 148
 149 /*! \class BamTools::BamAlignment
 150     \brief The main BAM alignment data structure.
 151
 152     Provides methods to query/modify BAM alignment data fields.
 153 */
 154 /*! \var BamAlignment::Name
 155     \brief read name
 156 */
 157 /*! \var BamAlignment::Length
 158     \brief length of query sequence
 159 */
 160 /*! \var BamAlignment::QueryBases
 161     \brief 'original' sequence (as reported from sequencing machine)
 162 */
 163 /*! \var BamAlignment::AlignedBases
 164     \brief 'aligned' sequence (includes any indels, padding, clipping)
 165 */
 166 /*! \var BamAlignment::Qualities
 167     \brief FASTQ qualities (ASCII characters, not numeric values)
 168 */
 169 /*! \var BamAlignment::TagData
 170     \brief tag data (use the provided methods to query/modify)
 171 */
 172 /*! \var BamAlignment::RefID
 173     \brief ID number for reference sequence
 174 */
 175 /*! \var BamAlignment::Position
 176     \brief position (0-based) where alignment starts
 177 */
 178 /*! \var BamAlignment::Bin
 179     \brief BAM (standard) index bin number for this alignment
 180 */
 181 /*! \var BamAlignment::MapQuality
 182     \brief mapping quality score
 183 */
 184 /*! \var BamAlignment::AlignmentFlag
 185     \brief alignment bit-flag (use the provided methods to query/modify)
 186 */
 187 /*! \var BamAlignment::CigarData
 188     \brief CIGAR operations for this alignment
 189 */
 190 /*! \var BamAlignment::MateRefID
 191     \brief ID number for reference sequence where alignment's mate was aligned
 192 */
 193 /*! \var BamAlignment::MatePosition
 194     \brief position (0-based) where alignment's mate starts
 195 */
 196 /*! \var BamAlignment::InsertSize
 197     \brief mate-pair insert size
 198 */
 199 /*! \var BamAlignment::Filename
 200     \brief name of BAM file which this alignment comes from
 201 */
 202
 203 /*! \fn BamAlignment::BamAlignment(void)
 204     \brief constructor
 205 */
 206 BamAlignment::BamAlignment(void)
 207     : RefID(-1)
 208     , Position(-1)
 209     , MateRefID(-1)
 210     , MatePosition(-1)
 211     , InsertSize(0)
 212 { }
 213
 214 /*! \fn BamAlignment::BamAlignment(const BamAlignment& other)
 215     \brief copy constructor
 216 */
 217 BamAlignment::BamAlignment(const BamAlignment& other)
 218     : Name(other.Name)
 219     , Length(other.Length)
 220     , QueryBases(other.QueryBases)
 221     , AlignedBases(other.AlignedBases)
 222     , Qualities(other.Qualities)
 223     , TagData(other.TagData)
 224     , RefID(other.RefID)
 225     , Position(other.Position)
 226     , Bin(other.Bin)
 227     , MapQuality(other.MapQuality)
 228     , AlignmentFlag(other.AlignmentFlag)
 229     , CigarData(other.CigarData)
 230     , MateRefID(other.MateRefID)
 231     , MatePosition(other.MatePosition)
 232     , InsertSize(other.InsertSize)
 233     , Filename(other.Filename)
 234     , SupportData(other.SupportData)
 235 { }
 236
 237 /*! \fn BamAlignment::~BamAlignment(void)
 238     \brief destructor
 239 */
 240 BamAlignment::~BamAlignment(void) { }
 241
 242 /*! \fn bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value)
 243     \brief Adds a field with string data to the BAM tags.
 244
 245     Does NOT modify an existing tag - use \link BamAlignment::EditTag() \endlink instead.
 246
 247     \param tag   2-character tag name
 248     \param type  1-character tag type (must be "Z" or "H")
 249     \param value string data to store
 250
 251     \return \c true if the \b new tag was added successfully
 252
 253     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 254         for more details on reserved tag names, supported tag types, etc.
 255 */
 256 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
 257
 258     // skip if core data not parsed
 259     if ( SupportData.HasCoreOnly ) return false;
 260
 261     // validate tag/type size & that type is OK for string value
 262     if ( !Internal::IsValidSize(tag, type) ) return false;
 263     if ( type.at(0) != Constants::BAM_TAG_TYPE_STRING &&
 264          type.at(0) != Constants::BAM_TAG_TYPE_HEX )
 265         return false;
 266
 267     // localize the tag data
 268     char* pTagData = (char*)TagData.data();
 269     const unsigned int tagDataLength = TagData.size();
 270     unsigned int numBytesParsed = 0;
 271
 272     // if tag already exists, return false
 273     // use EditTag explicitly instead
 274     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 275         return false;
 276
 277     // otherwise, copy tag data to temp buffer
 278     string newTag = tag + type + value;
 279     const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
 280     char originalTagData[newTagDataLength];
 281     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 282
 283     // append newTag
 284     strcat(originalTagData + tagDataLength, newTag.data());  // removes original null-term, appends newTag + null-term
 285
 286     // store temp buffer back in TagData
 287     const char* newTagData = (const char*)originalTagData;
 288     TagData.assign(newTagData, newTagDataLength);
 289
 290     // return success
 291     return true;
 292 }
 293
 294 /*! \fn bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value)
 295     \brief Adds a field with unsigned integer data to the BAM tags.
 296
 297     Does NOT modify an existing tag - use \link BamAlignment::EditTag() \endlink instead.
 298
 299     \param tag   2-character tag name
 300     \param type  1-character tag type (must NOT be "f", "Z", or "H")
 301     \param value unsigned int data to store
 302
 303     \return \c true if the \b new tag was added successfully
 304     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 305         for more details on reserved tag names, supported tag types, etc.
 306 */
 307 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
 308
 309     // skip if core data not parsed
 310     if ( SupportData.HasCoreOnly ) return false;
 311
 312     // validate tag/type size & that type is OK for uint32_t value
 313     if ( !Internal::IsValidSize(tag, type) ) return false;
 314     if ( type.at(0) == Constants::BAM_TAG_TYPE_FLOAT ||
 315          type.at(0) == Constants::BAM_TAG_TYPE_STRING ||
 316          type.at(0) == Constants::BAM_TAG_TYPE_HEX )
 317         return false;
 318
 319     // localize the tag data
 320     char* pTagData = (char*)TagData.data();
 321     const unsigned int tagDataLength = TagData.size();
 322     unsigned int numBytesParsed = 0;
 323
 324     // if tag already exists, return false
 325     // use EditTag explicitly instead
 326     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 327         return false;
 328
 329     // otherwise, convert value to string
 330     union { uint32_t value; char valueBuffer[sizeof(uint32_t)]; } un;
 331     un.value = value;
 332
 333     // copy original tag data to temp buffer
 334     string newTag = tag + type;
 335     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
 336     char originalTagData[newTagDataLength];
 337     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 338
 339     // append newTag
 340     strcat(originalTagData + tagDataLength, newTag.data());
 341     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(uint32_t));
 342
 343     // store temp buffer back in TagData
 344     const char* newTagData = (const char*)originalTagData;
 345     TagData.assign(newTagData, newTagDataLength);
 346
 347     // return success
 348     return true;
 349 }
 350
 351 /*! \fn bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value)
 352     \brief Adds a field with signed integer data to the BAM tags.
 353
 354     Does NOT modify an existing tag - use \link BamAlignment::EditTag() \endlink instead.
 355
 356     \param tag   2-character tag name
 357     \param type  1-character tag type (must NOT be "f", "Z", or "H")
 358     \param value signed int data to store
 359
 360     \return \c true if the \b new tag was added successfully
 361
 362     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 363         for more details on reserved tag names, supported tag types, etc.
 364 */
 365 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
 366     return AddTag(tag, type, (const uint32_t&)value);
 367 }
 368
 369 /*! \fn bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value)
 370     \brief Adds a field with floating-point data to the BAM tags.
 371
 372     Does NOT modify an existing tag - use \link BamAlignment::EditTag() \endlink instead.
 373
 374     \param tag  2-character tag name
 375     \param type  1-character tag type (must NOT be "Z" or "H")
 376     \param value float data to store
 377
 378     \return \c true if the \b new tag was added successfully
 379
 380     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 381         for more details on reserved tag names, supported tag types, etc.
 382 */
 383 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
 384
 385     // skip if core data not parsed
 386     if ( SupportData.HasCoreOnly ) return false;
 387
 388     // validate tag/type size & that type is OK for float value
 389     if ( !Internal::IsValidSize(tag, type) ) return false;
 390     if ( type.at(0) == Constants::BAM_TAG_TYPE_STRING ||
 391          type.at(0) == Constants::BAM_TAG_TYPE_HEX )
 392         return false;
 393
 394     // localize the tag data
 395     char* pTagData = (char*)TagData.data();
 396     const unsigned int tagDataLength = TagData.size();
 397     unsigned int numBytesParsed = 0;
 398
 399     // if tag already exists, return false
 400     // use EditTag explicitly instead
 401     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 402         return false;
 403
 404     // otherwise, convert value to string
 405     union { float value; char valueBuffer[sizeof(float)]; } un;
 406     un.value = value;
 407
 408     // copy original tag data to temp buffer
 409     string newTag = tag + type;
 410     const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
 411     char originalTagData[newTagDataLength];
 412     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 413
 414     // append newTag
 415     strcat(originalTagData + tagDataLength, newTag.data());
 416     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
 417
 418     // store temp buffer back in TagData
 419     const char* newTagData = (const char*)originalTagData;
 420     TagData.assign(newTagData, newTagDataLength);
 421
 422     // return success
 423     return true;
 424 }
 425
 426 /*! \fn bool BamAlignment::BuildCharData(void)
 427     \brief Populates alignment string fields (read name, bases, qualities, tag data).
 428
 429     An alignment retrieved using BamReader::GetNextAlignmentCore() lacks this data.
 430     Using that method makes parsing much quicker when only positional data is required.
 431
 432     However, if you later want to access the character data fields from such an alignment,
 433     use this method to populate those fields. Provides ability to do 'lazy evaluation' of
 434     alignment parsing.
 435
 436     \return \c true if character data populated successfully (or was already available to begin with)
 437 */
 438 bool BamAlignment::BuildCharData(void) {
 439
 440     // skip if char data already parsed
 441     if ( !SupportData.HasCoreOnly )
 442         return true;
 443
 444     // check system endianness
 445     bool IsBigEndian = BamTools::SystemIsBigEndian();
 446
 447     // calculate character lengths/offsets
 448     const unsigned int dataLength     = SupportData.BlockLength - Constants::BAM_CORE_SIZE;
 449     const unsigned int seqDataOffset  = SupportData.QueryNameLength + (SupportData.NumCigarOperations * 4);
 450     const unsigned int qualDataOffset = seqDataOffset + (SupportData.QuerySequenceLength+1)/2;
 451     const unsigned int tagDataOffset  = qualDataOffset + SupportData.QuerySequenceLength;
 452     const unsigned int tagDataLength  = dataLength - tagDataOffset;
 453
 454     // check offsets to see what char data exists
 455     const bool hasSeqData  = ( seqDataOffset  < dataLength );
 456     const bool hasQualData = ( qualDataOffset < dataLength );
 457     const bool hasTagData  = ( tagDataOffset  < dataLength );
 458
 459     // set up char buffers
 460     const char* allCharData = SupportData.AllCharData.data();
 461     const char* seqData     = ( hasSeqData  ? (((const char*)allCharData) + seqDataOffset)  : (const char*)0 );
 462     const char* qualData    = ( hasQualData ? (((const char*)allCharData) + qualDataOffset) : (const char*)0 );
 463           char* tagData     = ( hasTagData  ? (((char*)allCharData) + tagDataOffset)        : (char*)0 );
 464
 465     // store alignment name (relies on null char in name as terminator)
 466     Name.assign((const char*)(allCharData));
 467
 468     // save query sequence
 469     QueryBases.clear();
 470     if ( hasSeqData ) {
 471         QueryBases.reserve(SupportData.QuerySequenceLength);
 472         for (unsigned int i = 0; i < SupportData.QuerySequenceLength; ++i) {
 473             char singleBase = Constants::BAM_DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
 474             QueryBases.append(1, singleBase);
 475         }
 476     }
 477
 478     // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
 479     Qualities.clear();
 480     if ( hasQualData ) {
 481         Qualities.reserve(SupportData.QuerySequenceLength);
 482         for (unsigned int i = 0; i < SupportData.QuerySequenceLength; ++i) {
 483             char singleQuality = (char)(qualData[i]+33);
 484             Qualities.append(1, singleQuality);
 485         }
 486     }
 487
 488     // clear previous AlignedBases
 489     AlignedBases.clear();
 490
 491     // if QueryBases has data, build AlignedBases using CIGAR data
 492     // otherwise, AlignedBases will remain empty (this case IS allowed)
 493     if ( !QueryBases.empty() ) {
 494
 495         // resize AlignedBases
 496         AlignedBases.reserve(SupportData.QuerySequenceLength);
 497
 498         // iterate over CigarOps
 499         int k = 0;
 500         vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 501         vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 502         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 503             const CigarOp& op = (*cigarIter);
 504
 505             switch (op.Type) {
 506
 507                 // for 'M', 'I' - write bases
 508                 case (Constants::BAM_CIGAR_MATCH_CHAR) :
 509                 case (Constants::BAM_CIGAR_INS_CHAR) :
 510                     AlignedBases.append(QueryBases.substr(k, op.Length));
 511                     // fall through
 512
 513                 // for 'S' - soft clip, do not write bases
 514                 // but increment placeholder 'k'
 515                 case (Constants::BAM_CIGAR_SOFTCLIP_CHAR) :
 516                     k += op.Length;
 517                     break;
 518
 519                 // for 'D' - write gap character
 520                 case (Constants::BAM_CIGAR_DEL_CHAR) :
 521                     AlignedBases.append(op.Length, Constants::BAM_DNA_DEL);
 522                     break;
 523
 524                 // for 'P' - write padding character
 525                 case (Constants::BAM_CIGAR_PAD_CHAR) :
 526                     AlignedBases.append( op.Length, Constants::BAM_DNA_PAD );
 527                     break;
 528
 529                 // for 'N' - write N's, skip bases in original query sequence
 530                 case (Constants::BAM_CIGAR_REFSKIP_CHAR) :
 531                     AlignedBases.append( op.Length, Constants::BAM_DNA_N );
 532                     break;
 533
 534                 // for 'H' - hard clip, do nothing to AlignedBases, move to next op
 535                 case (Constants::BAM_CIGAR_HARDCLIP_CHAR) :
 536                     break;
 537
 538                 // shouldn't get here
 539                 default:
 540                     fprintf(stderr, "BamAlignment ERROR: invalid CIGAR operation type: %c\n", op.Type);
 541                     exit(1);
 542             }
 543         }
 544     }
 545
 546     // save tag data
 547     TagData.clear();
 548     if ( hasTagData ) {
 549         if ( IsBigEndian ) {
 550             int i = 0;
 551             while ( (unsigned int)i < tagDataLength ) {
 552
 553                 i += Constants::BAM_TAG_TAGSIZE;  // skip tag chars (e.g. "RG", "NM", etc.)
 554                 const char type = tagData[i];     // get tag type at position i
 555                 ++i;                              // move i past tag type
 556
 557                 switch (type) {
 558
 559                     case(Constants::BAM_TAG_TYPE_ASCII) :
 560                     case(Constants::BAM_TAG_TYPE_INT8)  :
 561                     case(Constants::BAM_TAG_TYPE_UINT8) :
 562                         ++i;
 563                         break;
 564
 565                     case(Constants::BAM_TAG_TYPE_INT16)  :
 566                     case(Constants::BAM_TAG_TYPE_UINT16) :
 567                         BamTools::SwapEndian_16p(&tagData[i]);
 568                         i += sizeof(uint16_t);
 569                         break;
 570
 571                     case(Constants::BAM_TAG_TYPE_FLOAT)  :
 572                     case(Constants::BAM_TAG_TYPE_INT32)  :
 573                     case(Constants::BAM_TAG_TYPE_UINT32) :
 574                         BamTools::SwapEndian_32p(&tagData[i]);
 575                         i += sizeof(uint32_t);
 576                         break;
 577
 578                     case(Constants::BAM_TAG_TYPE_HEX) :
 579                     case(Constants::BAM_TAG_TYPE_STRING) :
 580                         // no endian swapping necessary for hex-string/string data
 581                         while (tagData[i]) { ++i; }
 582                         // increment one more for null terminator
 583                         ++i;
 584                         break;
 585
 586                     // shouldn't get here
 587                     default :
 588                         fprintf(stderr, "BamAlignment ERROR: invalid tag value type: %c\n", type);
 589                         exit(1);
 590                 }
 591             }
 592         }
 593
 594         // store tagData in alignment
 595         TagData.resize(tagDataLength);
 596         memcpy((char*)TagData.data(), tagData, tagDataLength);
 597     }
 598
 599     // clear the core-only flag
 600     SupportData.HasCoreOnly = false;
 601
 602     // return success
 603     return true;
 604 }
 605
 606 /*! \fn bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value)
 607     \brief Edits a BAM tag field containing string data.
 608
 609     If \a tag does not exist, a new entry is created.
 610
 611     \param tag   2-character tag name
 612     \param type  1-character tag type (must be "Z" or "H")
 613     \param value string data to store
 614
 615     \return \c true if the tag was modified/created successfully
 616
 617     \sa BamAlignment::RemoveTag()
 618     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 619         for more details on reserved tag names, supported tag types, etc.
 620 */
 621 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
 622
 623     // skip if core data not parsed
 624     if ( SupportData.HasCoreOnly ) return false;
 625
 626     // validate tag/type size & that type is OK for string value
 627     if ( !Internal::IsValidSize(tag, type) ) return false;
 628     if ( type.at(0) != Constants::BAM_TAG_TYPE_STRING &&
 629          type.at(0) != Constants::BAM_TAG_TYPE_HEX )
 630         return false;
 631
 632     // localize the tag data
 633     char* pOriginalTagData = (char*)TagData.data();
 634     char* pTagData = pOriginalTagData;
 635     const unsigned int originalTagDataLength = TagData.size();
 636
 637     unsigned int newTagDataLength = 0;
 638     unsigned int numBytesParsed = 0;
 639
 640     // if tag found, store data in readGroup, return success
 641     if ( Internal::FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 642
 643         // make sure array is more than big enough
 644         char newTagData[originalTagDataLength + value.size()];
 645
 646         // copy original tag data up til desired tag
 647         const unsigned int beginningTagDataLength = numBytesParsed;
 648         newTagDataLength += beginningTagDataLength;
 649         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 650
 651         // copy new VALUE in place of current tag data
 652         const unsigned int dataLength = strlen(value.c_str());
 653         memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
 654
 655         // skip to next tag (if tag for removal is last, return true)
 656         const char* pTagStorageType = pTagData - 1;
 657         if ( !Internal::SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) )
 658             return true;
 659
 660         // copy everything from current tag (the next one after tag for removal) to end
 661         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 662         const unsigned int endTagOffset      = beginningTagDataLength + dataLength + 1;
 663         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 664         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 665
 666         // ensure null-terminator
 667         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 668
 669         // save new tag data
 670         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 671         return true;
 672     }
 673
 674     // tag not found, attempt AddTag
 675     else return AddTag(tag, type, value);
 676 }
 677
 678 /*! \fn bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value)
 679     \brief Edits a BAM tag field containing unsigned integer data.
 680
 681     If \a tag does not exist, a new entry is created.
 682
 683     \param tag   2-character tag name
 684     \param type  1-character tag type (must NOT be "f", "Z", or "H")
 685     \param value unsigned integer data to store
 686
 687     \return \c true if the tag was modified/created successfully
 688
 689     \sa BamAlignment::RemoveTag()
 690     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 691         for more details on reserved tag names, supported tag types, etc.
 692 */
 693 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
 694
 695     // skip if core data not parsed
 696     if ( SupportData.HasCoreOnly ) return false;
 697
 698     // validate tag/type size & that type is OK for uint32_t value
 699     if ( !Internal::IsValidSize(tag, type) ) return false;
 700     if ( type.at(0) == Constants::BAM_TAG_TYPE_FLOAT ||
 701          type.at(0) == Constants::BAM_TAG_TYPE_STRING ||
 702          type.at(0) == Constants::BAM_TAG_TYPE_HEX )
 703         return false;
 704
 705     // localize the tag data
 706     char* pOriginalTagData = (char*)TagData.data();
 707     char* pTagData = pOriginalTagData;
 708     const unsigned int originalTagDataLength = TagData.size();
 709
 710     unsigned int newTagDataLength = 0;
 711     unsigned int numBytesParsed = 0;
 712
 713     // if tag found, store data in readGroup, return success
 714     if ( Internal::FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 715
 716         // make sure array is more than big enough
 717         char newTagData[originalTagDataLength + sizeof(value)];
 718
 719         // copy original tag data up til desired tag
 720         const unsigned int beginningTagDataLength = numBytesParsed;
 721         newTagDataLength += beginningTagDataLength;
 722         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 723
 724         // copy new VALUE in place of current tag data
 725         union { uint32_t value; char valueBuffer[sizeof(uint32_t)]; } un;
 726         un.value = value;
 727         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(uint32_t));
 728
 729         // skip to next tag (if tag for removal is last, return true)
 730         const char* pTagStorageType = pTagData - 1;
 731         if ( !Internal::SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) )
 732             return true;
 733
 734         // copy everything from current tag (the next one after tag for removal) to end
 735         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 736         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(uint32_t);
 737         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 738         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 739
 740         // ensure null-terminator
 741         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 742
 743         // save new tag data
 744         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 745         return true;
 746     }
 747
 748     // tag not found, attempt AddTag
 749     else return AddTag(tag, type, value);
 750 }
 751
 752 /*! \fn bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value)
 753     \brief Edits a BAM tag field containing signed integer data.
 754
 755     If \a tag does not exist, a new entry is created.
 756
 757     \param tag   2-character tag name
 758     \param type  1-character tag type (must NOT be "f", "Z", or "H")
 759     \param value signed integer data to store
 760
 761     \return \c true if the tag was modified/created successfully
 762
 763     \sa BamAlignment::RemoveTag()
 764     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 765         for more details on reserved tag names, supported tag types, etc.
 766 */
 767 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
 768     return EditTag(tag, type, (const uint32_t&)value);
 769 }
 770
 771 /*! \fn bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value)
 772     \brief Edits a BAM tag field containing floating-point data.
 773
 774     If \a tag does not exist, a new entry is created.
 775
 776     \param tag   2-character tag name
 777     \param type  1-character tag type (must NOT be "Z" or "H")
 778     \param value float data to store
 779
 780     \return \c true if the tag was modified/created successfully
 781
 782     \sa BamAlignment::RemoveTag()
 783     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
 784         for more details on reserved tag names, supported tag types, etc.
 785 */
 786 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
 787
 788     // skip if core data not parsed
 789     if ( SupportData.HasCoreOnly ) return false;
 790
 791     // validate tag/type size & that type is OK for float value
 792     if ( !Internal::IsValidSize(tag, type) ) return false;
 793     if ( type.at(0) == Constants::BAM_TAG_TYPE_STRING ||
 794          type.at(0) == Constants::BAM_TAG_TYPE_HEX )
 795         return false;
 796
 797      // localize the tag data
 798     char* pOriginalTagData = (char*)TagData.data();
 799     char* pTagData = pOriginalTagData;
 800     const unsigned int originalTagDataLength = TagData.size();
 801
 802     unsigned int newTagDataLength = 0;
 803     unsigned int numBytesParsed = 0;
 804
 805     // if tag found, store data in readGroup, return success
 806     if ( Internal::FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
 807
 808         // make sure array is more than big enough
 809         char newTagData[originalTagDataLength + sizeof(value)];
 810
 811         // copy original tag data up til desired tag
 812         const unsigned int beginningTagDataLength = numBytesParsed;
 813         newTagDataLength += beginningTagDataLength;
 814         memcpy(newTagData, pOriginalTagData, numBytesParsed);
 815
 816         // copy new VALUE in place of current tag data
 817         union { float value; char valueBuffer[sizeof(float)]; } un;
 818         un.value = value;
 819         memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
 820
 821         // skip to next tag (if tag for removal is last, return true)
 822         const char* pTagStorageType = pTagData - 1;
 823         if ( !Internal::SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) )
 824             return true;
 825
 826         // copy everything from current tag (the next one after tag for removal) to end
 827         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 828         const unsigned int endTagOffset      = beginningTagDataLength + sizeof(float);
 829         const unsigned int endTagDataLength  = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 830         memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
 831
 832         // ensure null-terminator
 833         newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
 834
 835         // save new tag data
 836         TagData.assign(newTagData, endTagOffset + endTagDataLength);
 837         return true;
 838     }
 839
 840     // tag not found, attempt AddTag
 841     else return AddTag(tag, type, value);
 842 }
 843
 844 /*! \fn bool BamAlignment::GetEditDistance(uint32_t& editDistance) const
 845     \brief Retrieves value of edit distance tag ("NM").
 846
 847     \deprecated Instead use BamAlignment::GetTag()
 848         \code
 849             BamAlignment::GetTag("NM", editDistance);
 850         \endcode
 851
 852     \param editDistance destination for retrieved value
 853
 854     \return \c true if found
 855 */
 856 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
 857     return GetTag("NM", (uint32_t&)editDistance);
 858 }
 859
 860 /*! \fn int BamAlignment::GetEndPosition(bool usePadded = false, bool zeroBased = true) const
 861     \brief Calculates alignment end position, based on starting position and CIGAR data.
 862
 863     \param usePadded Inserted bases affect reported position. Default is false, so that reported
 864                      position stays 'sync-ed' with reference coordinates.
 865     \param zeroBased Return (BAM standard) 0-based coordinate. Setting this to false can be useful
 866                      when using BAM data with half-open formats (e.g. BED).
 867
 868     \return alignment end position
 869 */
 870 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {
 871
 872     // initialize alignment end to starting position
 873     int alignEnd = Position;
 874
 875     // iterate over cigar operations
 876     vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 877     vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 878     for ( ; cigarIter != cigarEnd; ++cigarIter) {
 879         const char cigarType = (*cigarIter).Type;
 880         const uint32_t& cigarLength = (*cigarIter).Length;
 881
 882         if ( cigarType == Constants::BAM_CIGAR_MATCH_CHAR ||
 883              cigarType == Constants::BAM_CIGAR_DEL_CHAR ||
 884              cigarType == Constants::BAM_CIGAR_REFSKIP_CHAR )
 885             alignEnd += cigarLength;
 886         else if ( usePadded && cigarType == Constants::BAM_CIGAR_INS_CHAR )
 887             alignEnd += cigarLength;
 888     }
 889
 890     // adjust for zero-based coordinates, if requested
 891     if ( zeroBased ) alignEnd -= 1;
 892
 893     // return result
 894     return alignEnd;
 895 }
 896
 897 /*! \fn bool BamAlignment::GetReadGroup(std::string& readGroup) const
 898     \brief Retrieves value of read group tag ("RG").
 899
 900     \deprecated Instead use BamAlignment::GetTag()
 901         \code
 902             BamAlignment::GetTag("RG", readGroup);
 903         \endcode
 904
 905     \param readGroup destination for retrieved value
 906
 907     \return \c true if found
 908 */
 909 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
 910     return GetTag("RG", readGroup);
 911 }
 912
 913 /*! \fn bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const
 914     \brief Retrieves the string value associated with a BAM tag.
 915
 916     \param tag         2-character tag name
 917     \param destination destination for retrieved value
 918
 919     \return \c true if found
 920 */
 921 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
 922
 923     // make sure tag data exists
 924     if ( SupportData.HasCoreOnly || TagData.empty() )
 925         return false;
 926
 927     // localize the tag data
 928     char* pTagData = (char*)TagData.data();
 929     const unsigned int tagDataLength = TagData.size();
 930     unsigned int numBytesParsed = 0;
 931
 932     // if tag found, store data in readGroup, return success
 933     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 934         const unsigned int dataLength = strlen(pTagData);
 935         destination.clear();
 936         destination.resize(dataLength);
 937         memcpy( (char*)destination.data(), pTagData, dataLength );
 938         return true;
 939     }
 940
 941     // tag not found, return failure
 942     return false;
 943 }
 944
 945 /*! \fn bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const
 946     \brief Retrieves the unsigned integer value associated with a BAM tag.
 947
 948     \param tag         2-character tag name
 949     \param destination destination for retrieved value
 950
 951     \return \c true if found
 952 */
 953 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
 954
 955     // make sure tag data exists
 956     if ( SupportData.HasCoreOnly || TagData.empty() )
 957         return false;
 958
 959     // localize the tag data
 960     char* pTagData = (char*)TagData.data();
 961     const unsigned int tagDataLength = TagData.size();
 962     unsigned int numBytesParsed = 0;
 963
 964     // if tag found, determine data byte-length, store data in readGroup, return success
 965     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
 966
 967         // determine data byte-length
 968         const char type = *(pTagData - 1);
 969         int destinationLength = 0;
 970         switch (type) {
 971
 972             // 1 byte data
 973             case (Constants::BAM_TAG_TYPE_ASCII) :
 974             case (Constants::BAM_TAG_TYPE_INT8)  :
 975             case (Constants::BAM_TAG_TYPE_UINT8) :
 976                 destinationLength = 1;
 977                 break;
 978
 979             // 2 byte data
 980             case (Constants::BAM_TAG_TYPE_INT16)  :
 981             case (Constants::BAM_TAG_TYPE_UINT16) :
 982                 destinationLength = 2;
 983                 break;
 984
 985             // 4 byte data
 986             case (Constants::BAM_TAG_TYPE_INT32)  :
 987             case (Constants::BAM_TAG_TYPE_UINT32) :
 988                 destinationLength = 4;
 989                 break;
 990
 991             // unsupported type for integer destination (float or var-length strings)
 992             case (Constants::BAM_TAG_TYPE_FLOAT)  :
 993             case (Constants::BAM_TAG_TYPE_STRING) :
 994             case (Constants::BAM_TAG_TYPE_HEX)    :
 995                 fprintf(stderr, "BamAlignment ERROR: cannot store tag of type %c in integer destination\n", type);
 996                 return false;
 997
 998             // unknown tag type
 999             default:
1000                 fprintf(stderr, "BamAlignment ERROR: unknown tag type encountered: [%c]\n", type);
1001                 return false;
1002         }
1003
1004         // store in destination
1005         destination = 0;
1006         memcpy(&destination, pTagData, destinationLength);
1007         return true;
1008     }
1009
1010     // tag not found, return failure
1011     return false;
1012 }
1013
1014 /*! \fn bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const
1015     \brief Retrieves the signed integer value associated with a BAM tag.
1016
1017     \param tag         2-character tag name
1018     \param destination destination for retrieved value
1019
1020     \return \c true if found
1021 */
1022 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
1023     return GetTag(tag, (uint32_t&)destination);
1024 }
1025
1026 /*! \fn bool BamAlignment::GetTag(const std::string& tag, float& destination) const
1027     \brief Retrieves the floating-point value associated with a BAM tag.
1028
1029     \param tag         2-character tag name
1030     \param destination destination for retrieved value
1031
1032     \return \c true if found
1033 */
1034 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
1035
1036     // make sure tag data exists
1037     if ( SupportData.HasCoreOnly || TagData.empty() )
1038         return false;
1039
1040     // localize the tag data
1041     char* pTagData = (char*)TagData.data();
1042     const unsigned int tagDataLength = TagData.size();
1043     unsigned int numBytesParsed = 0;
1044
1045     // if tag found, determine data byte-length, store data in readGroup, return success
1046     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
1047
1048         // determine data byte-length
1049         const char type = *(pTagData - 1);
1050         int destinationLength = 0;
1051         switch (type) {
1052
1053             // 1 byte data
1054             case (Constants::BAM_TAG_TYPE_ASCII) :
1055             case (Constants::BAM_TAG_TYPE_INT8)  :
1056             case (Constants::BAM_TAG_TYPE_UINT8) :
1057                 destinationLength = 1;
1058                 break;
1059
1060             // 2 byte data
1061             case (Constants::BAM_TAG_TYPE_INT16)  :
1062             case (Constants::BAM_TAG_TYPE_UINT16) :
1063                 destinationLength = 2;
1064                 break;
1065
1066             // 4 byte data
1067             case (Constants::BAM_TAG_TYPE_FLOAT)  :
1068             case (Constants::BAM_TAG_TYPE_INT32)  :
1069             case (Constants::BAM_TAG_TYPE_UINT32) :
1070                 destinationLength = 4;
1071                 break;
1072
1073             // unsupported type (var-length strings)
1074             case (Constants::BAM_TAG_TYPE_STRING) :
1075             case (Constants::BAM_TAG_TYPE_HEX)    :
1076                 fprintf(stderr, "BamAlignment ERROR: cannot store tag of type %c in float destination\n", type);
1077                 return false;
1078
1079             // unknown tag type
1080             default:
1081                 fprintf(stderr, "BamAlignment ERROR: unknown tag type encountered: [%c]\n", type);
1082                 return false;
1083         }
1084
1085         // store in destination
1086         destination = 0.0;
1087         memcpy(&destination, pTagData, destinationLength);
1088         return true;
1089     }
1090
1091     // tag not found, return failure
1092     return false;
1093 }
1094
1095 /*! \fn bool BamAlignment::GetTagType(const std::string& tag, char& type) const
1096     \brief Retrieves the BAM tag type-code associated with requested tag name.
1097
1098     \param tag  2-character tag name
1099     \param type destination for the retrieved (1-character) tag type
1100
1101     \return \c true if found
1102
1103     \sa http://samtools.sourceforge.net/SAM-1.3.pdf
1104         for more details on reserved tag names, supported tag types, etc.
1105 */
1106 bool BamAlignment::GetTagType(const std::string& tag, char& type) const {
1107
1108     // make sure tag data exists
1109     if ( SupportData.HasCoreOnly || TagData.empty() )
1110         return false;
1111
1112     // localize the tag data
1113     char* pTagData = (char*)TagData.data();
1114     const unsigned int tagDataLength = TagData.size();
1115     unsigned int numBytesParsed = 0;
1116
1117     // lookup tag
1118     if ( Internal::FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
1119
1120         // retrieve tag type code
1121         type = *(pTagData - 1);
1122
1123         // validate that type is a proper BAM tag type
1124         switch (type) {
1125             case (Constants::BAM_TAG_TYPE_ASCII)  :
1126             case (Constants::BAM_TAG_TYPE_INT8)   :
1127             case (Constants::BAM_TAG_TYPE_UINT8)  :
1128             case (Constants::BAM_TAG_TYPE_INT16)  :
1129             case (Constants::BAM_TAG_TYPE_UINT16) :
1130             case (Constants::BAM_TAG_TYPE_INT32)  :
1131             case (Constants::BAM_TAG_TYPE_UINT32) :
1132             case (Constants::BAM_TAG_TYPE_FLOAT)  :
1133             case (Constants::BAM_TAG_TYPE_STRING) :
1134             case (Constants::BAM_TAG_TYPE_HEX)    :
1135                 return true;
1136
1137             // unknown tag type
1138             default:
1139                 fprintf(stderr, "BamAlignment ERROR: unknown tag type encountered: [%c]\n", type);
1140                 return false;
1141         }
1142     }
1143
1144     // tag not found, return failure
1145     return false;
1146 }
1147
1148 /*! \fn bool BamAlignment::IsDuplicate(void) const
1149     \return \c true if this read is a PCR duplicate
1150 */
1151 bool BamAlignment::IsDuplicate(void) const {
1152     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_DUPLICATE) != 0 );
1153 }
1154
1155 /*! \fn bool BamAlignment::IsFailedQC(void) const
1156     \return \c true if this read failed quality control
1157 */
1158 bool BamAlignment::IsFailedQC(void) const {
1159     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_QC_FAILED) != 0 );
1160 }
1161
1162 /*! \fn bool BamAlignment::IsFirstMate(void) const
1163     \return \c true if alignment is first mate on paired-end read
1164 */
1165 bool BamAlignment::IsFirstMate(void) const {
1166     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_1) != 0 );
1167 }
1168
1169 /*! \fn bool BamAlignment::IsMapped(void) const
1170     \return \c true if alignment is mapped
1171 */
1172 bool BamAlignment::IsMapped(void) const {
1173     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_UNMAPPED) == 0 );
1174 }
1175
1176 /*! \fn bool BamAlignment::IsMateMapped(void) const
1177     \return \c true if alignment's mate is mapped
1178 */
1179 bool BamAlignment::IsMateMapped(void) const {
1180     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_UNMAPPED) == 0 );
1181 }
1182
1183 /*! \fn bool BamAlignment::IsMateReverseStrand(void) const
1184     \return \c true if alignment's mate mapped to reverse strand
1185 */
1186 bool BamAlignment::IsMateReverseStrand(void) const {
1187     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND) != 0 );
1188 }
1189
1190 /*! \fn bool BamAlignment::IsPaired(void) const
1191     \return \c true if alignment part of paired-end read
1192 */
1193 bool BamAlignment::IsPaired(void) const {
1194     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PAIRED) != 0 );
1195 }
1196
1197 /*! \fn bool BamAlignment::IsPrimaryAlignment(void) const
1198     \return \c true if reported position is primary alignment
1199 */
1200 bool BamAlignment::IsPrimaryAlignment(void) const  {
1201     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_SECONDARY) == 0 );
1202 }
1203
1204 /*! \fn bool BamAlignment::IsProperPair(void) const
1205     \return \c true if alignment is part of read that satisfied paired-end resolution
1206 */
1207 bool BamAlignment::IsProperPair(void) const {
1208     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PROPER_PAIR) != 0 );
1209 }
1210
1211 /*! \fn bool BamAlignment::IsReverseStrand(void) const
1212     \return \c true if alignment mapped to reverse strand
1213 */
1214 bool BamAlignment::IsReverseStrand(void) const {
1215     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_REVERSE_STRAND) != 0 );
1216 }
1217
1218 /*! \fn bool BamAlignment::IsSecondMate(void) const
1219     \return \c true if alignment is second mate on read
1220 */
1221 bool BamAlignment::IsSecondMate(void) const {
1222     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_2) != 0 );
1223 }
1224
1225 /*! \fn bool BamAlignment::RemoveTag(const std::string& tag)
1226     \brief Removes field from BAM tags.
1227
1228     \return \c true if tag was removed successfully (or didn't exist before)
1229 */
1230 bool BamAlignment::RemoveTag(const std::string& tag) {
1231
1232     // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
1233     // also, return false if no data present to remove
1234     if ( SupportData.HasCoreOnly || TagData.empty() )
1235         return false;
1236
1237     // localize the tag data
1238     char* pOriginalTagData = (char*)TagData.data();
1239     char* pTagData = pOriginalTagData;
1240     const unsigned int originalTagDataLength = TagData.size();
1241     unsigned int newTagDataLength = 0;
1242     unsigned int numBytesParsed = 0;
1243
1244     // if tag found, store data in readGroup, return success
1245     if ( Internal::FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
1246
1247         char newTagData[originalTagDataLength];
1248
1249         // copy original tag data up til desired tag
1250         pTagData       -= 3;
1251         numBytesParsed -= 3;
1252         const unsigned int beginningTagDataLength = numBytesParsed;
1253         newTagDataLength += beginningTagDataLength;
1254         memcpy(newTagData, pOriginalTagData, numBytesParsed);
1255
1256         // skip to next tag (if tag for removal is last, return true)
1257         const char* pTagStorageType = pTagData + 2;
1258         pTagData       += 3;
1259         numBytesParsed += 3;
1260         if ( !Internal::SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) )
1261             return true;
1262
1263         // copy everything from current tag (the next one after tag for removal) to end
1264         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
1265         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
1266         memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
1267
1268         // save new tag data
1269         TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
1270         return true;
1271     }
1272
1273     // tag not found, no removal - return failure
1274     return false;
1275 }
1276
1277 /*! \fn void BamAlignment::SetIsDuplicate(bool ok)
1278     \brief Sets value of "PCR duplicate" flag to \a ok.
1279 */
1280 void BamAlignment::SetIsDuplicate(bool ok) {
1281     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_DUPLICATE;
1282     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_DUPLICATE;
1283 }
1284
1285 /*! \fn void BamAlignment::SetIsFailedQC(bool ok)
1286     \brief Sets "failed quality control" flag to \a ok.
1287 */
1288 void BamAlignment::SetIsFailedQC(bool ok) {
1289     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_QC_FAILED;
1290     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_QC_FAILED;
1291 }
1292
1293 /*! \fn void BamAlignment::SetIsFirstMate(bool ok)
1294     \brief Sets "alignment is first mate" flag to \a ok.
1295 */
1296 void BamAlignment::SetIsFirstMate(bool ok) {
1297     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_1;
1298     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_1;
1299 }
1300
1301 /*! \fn void BamAlignment::SetIsMapped(bool ok)
1302     \brief Sets "alignment is mapped" flag to \a ok.
1303 */
1304 void BamAlignment::SetIsMapped(bool ok) {
1305     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_UNMAPPED;
1306     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_UNMAPPED;
1307 }
1308
1309 /*! \fn void BamAlignment::SetIsMateMapped(bool ok)
1310     \brief Sets "alignment's mate is mapped" flag to \a ok.
1311 */
1312 void BamAlignment::SetIsMateMapped(bool ok) {
1313     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
1314     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
1315 }
1316
1317 /*! \fn void BamAlignment::SetIsMateUnmapped(bool ok)
1318     \brief Complement of using SetIsMateMapped().
1319     \deprecated For sake of symmetry with the query methods
1320     \sa IsMateMapped(), SetIsMateMapped()
1321 */
1322 void BamAlignment::SetIsMateUnmapped(bool ok) {
1323     SetIsMateMapped(!ok);
1324 }
1325
1326 /*! \fn void BamAlignment::SetIsMateReverseStrand(bool ok)
1327     \brief Sets "alignment's mate mapped to reverse strand" flag to \a ok.
1328 */
1329 void BamAlignment::SetIsMateReverseStrand(bool ok) {
1330     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
1331     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
1332 }
1333
1334 /*! \fn void BamAlignment::SetIsPaired(bool ok)
1335     \brief Sets "alignment part of paired-end read" flag to \a ok.
1336 */
1337 void BamAlignment::SetIsPaired(bool ok) {
1338     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PAIRED;
1339     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PAIRED;
1340 }
1341
1342 /*! \fn void BamAlignment::SetIsPrimaryAlignment(bool ok)
1343     \brief Sets "position is primary alignment" flag to \a ok.
1344 */
1345 void BamAlignment::SetIsPrimaryAlignment(bool ok) {
1346     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_SECONDARY;
1347     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_SECONDARY;
1348 }
1349
1350 /*! \fn void BamAlignment::SetIsProperPair(bool ok)
1351     \brief Sets "alignment is part of read that satisfied paired-end resolution" flag to \a ok.
1352 */
1353 void BamAlignment::SetIsProperPair(bool ok) {
1354     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PROPER_PAIR;
1355     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PROPER_PAIR;
1356 }
1357
1358 /*! \fn void BamAlignment::SetIsReverseStrand(bool ok)
1359     \brief Sets "alignment mapped to reverse strand" flag to \a ok.
1360 */
1361 void BamAlignment::SetIsReverseStrand(bool ok) {
1362     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_REVERSE_STRAND;
1363     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_REVERSE_STRAND;
1364 }
1365
1366 /*! \fn void BamAlignment::SetIsSecondaryAlignment(bool ok)
1367     \brief Complement of using SetIsPrimaryAlignment().
1368     \deprecated For sake of symmetry with the query methods
1369     \sa IsPrimaryAlignment(), SetIsPrimaryAlignment()
1370 */
1371 void BamAlignment::SetIsSecondaryAlignment(bool ok) {
1372     SetIsPrimaryAlignment(!ok);
1373 }
1374
1375 /*! \fn void BamAlignment::SetIsSecondMate(bool ok)
1376     \brief Sets "alignment is second mate on read" flag to \a ok.
1377 */
1378 void BamAlignment::SetIsSecondMate(bool ok) {
1379     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_2;
1380     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_2;
1381 }
1382
1383 /*! \fn void BamAlignment::SetIsUnmapped(bool ok)
1384     \brief Complement of using SetIsMapped().
1385     \deprecated For sake of symmetry with the query methods
1386     \sa IsMapped(), SetIsMapped()
1387 */
1388 void BamAlignment::SetIsUnmapped(bool ok) {
1389     SetIsMapped(!ok);
1390 }