src/api/BamAlignment.cpp

   1 // ***************************************************************************
   2 // BamAlignment.cpp (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 22 February 2012 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides the BamAlignment data structure
   8 // ***************************************************************************
   9
  10 #include "api/BamAlignment.h"
  11 #include "api/BamConstants.h"
  12 using namespace BamTools;
  13 using namespace std;
  14
  15 /*! \class BamTools::BamAlignment
  16     \brief The main BAM alignment data structure.
  17
  18     Provides methods to query/modify BAM alignment data fields.
  19 */
  20 /*! \var BamAlignment::Name
  21     \brief read name
  22 */
  23 /*! \var BamAlignment::Length
  24     \brief length of query sequence
  25 */
  26 /*! \var BamAlignment::QueryBases
  27     \brief 'original' sequence (as reported from sequencing machine)
  28 */
  29 /*! \var BamAlignment::AlignedBases
  30     \brief 'aligned' sequence (includes any indels, padding, clipping)
  31 */
  32 /*! \var BamAlignment::Qualities
  33     \brief FASTQ qualities (ASCII characters, not numeric values)
  34 */
  35 /*! \var BamAlignment::TagData
  36     \brief tag data (use the provided methods to query/modify)
  37 */
  38 /*! \var BamAlignment::RefID
  39     \brief ID number for reference sequence
  40 */
  41 /*! \var BamAlignment::Position
  42     \brief position (0-based) where alignment starts
  43 */
  44 /*! \var BamAlignment::Bin
  45     \brief BAM (standard) index bin number for this alignment
  46 */
  47 /*! \var BamAlignment::MapQuality
  48     \brief mapping quality score
  49 */
  50 /*! \var BamAlignment::AlignmentFlag
  51     \brief alignment bit-flag (use the provided methods to query/modify)
  52 */
  53 /*! \var BamAlignment::CigarData
  54     \brief CIGAR operations for this alignment
  55 */
  56 /*! \var BamAlignment::MateRefID
  57     \brief ID number for reference sequence where alignment's mate was aligned
  58 */
  59 /*! \var BamAlignment::MatePosition
  60     \brief position (0-based) where alignment's mate starts
  61 */
  62 /*! \var BamAlignment::InsertSize
  63     \brief mate-pair insert size
  64 */
  65 /*! \var BamAlignment::Filename
  66     \brief name of BAM file which this alignment comes from
  67 */
  68
  69 /*! \fn BamAlignment::BamAlignment(void)
  70     \brief constructor
  71 */
  72 BamAlignment::BamAlignment(void)
  73     : Length(0)
  74     , RefID(-1)
  75     , Position(-1)
  76     , Bin(0)
  77     , MapQuality(0)
  78     , AlignmentFlag(0)
  79     , MateRefID(-1)
  80     , MatePosition(-1)
  81     , InsertSize(0)
  82 { }
  83
  84 /*! \fn BamAlignment::BamAlignment(const BamAlignment& other)
  85     \brief copy constructor
  86 */
  87 BamAlignment::BamAlignment(const BamAlignment& other)
  88     : Name(other.Name)
  89     , Length(other.Length)
  90     , QueryBases(other.QueryBases)
  91     , AlignedBases(other.AlignedBases)
  92     , Qualities(other.Qualities)
  93     , TagData(other.TagData)
  94     , RefID(other.RefID)
  95     , Position(other.Position)
  96     , Bin(other.Bin)
  97     , MapQuality(other.MapQuality)
  98     , AlignmentFlag(other.AlignmentFlag)
  99     , CigarData(other.CigarData)
 100     , MateRefID(other.MateRefID)
 101     , MatePosition(other.MatePosition)
 102     , InsertSize(other.InsertSize)
 103     , Filename(other.Filename)
 104     , SupportData(other.SupportData)
 105 { }
 106
 107 /*! \fn BamAlignment::~BamAlignment(void)
 108     \brief destructor
 109 */
 110 BamAlignment::~BamAlignment(void) { }
 111
 112 /*! \fn bool BamAlignment::BuildCharData(void)
 113     \brief Populates alignment string fields (read name, bases, qualities, tag data).
 114
 115     An alignment retrieved using BamReader::GetNextAlignmentCore() lacks this data.
 116     Using that method makes parsing much quicker when only positional data is required.
 117
 118     However, if you later want to access the character data fields from such an alignment,
 119     use this method to populate those fields. Provides ability to do 'lazy evaluation' of
 120     alignment parsing.
 121
 122     \return \c true if character data populated successfully (or was already available to begin with)
 123 */
 124 bool BamAlignment::BuildCharData(void) {
 125
 126     // skip if char data already parsed
 127     if ( !SupportData.HasCoreOnly )
 128         return true;
 129
 130     // check system endianness
 131     bool IsBigEndian = BamTools::SystemIsBigEndian();
 132
 133     // calculate character lengths/offsets
 134     const unsigned int dataLength     = SupportData.BlockLength - Constants::BAM_CORE_SIZE;
 135     const unsigned int seqDataOffset  = SupportData.QueryNameLength + (SupportData.NumCigarOperations*4);
 136     const unsigned int qualDataOffset = seqDataOffset + (SupportData.QuerySequenceLength+1)/2;
 137     const unsigned int tagDataOffset  = qualDataOffset + SupportData.QuerySequenceLength;
 138     const unsigned int tagDataLength  = dataLength - tagDataOffset;
 139
 140     // check offsets to see what char data exists
 141     const bool hasSeqData  = ( seqDataOffset  < dataLength );
 142     const bool hasQualData = ( qualDataOffset < dataLength );
 143     const bool hasTagData  = ( tagDataOffset  < dataLength );
 144
 145     // set up char buffers
 146     const char* allCharData = SupportData.AllCharData.data();
 147     const char* seqData     = ( hasSeqData  ? (((const char*)allCharData) + seqDataOffset)  : (const char*)0 );
 148     const char* qualData    = ( hasQualData ? (((const char*)allCharData) + qualDataOffset) : (const char*)0 );
 149           char* tagData     = ( hasTagData  ? (((char*)allCharData) + tagDataOffset)        : (char*)0 );
 150
 151     // store alignment name (relies on null char in name as terminator)
 152     Name.assign((const char*)(allCharData));
 153
 154     // save query sequence
 155     QueryBases.clear();
 156     if ( hasSeqData ) {
 157         QueryBases.reserve(SupportData.QuerySequenceLength);
 158         for ( size_t i = 0; i < SupportData.QuerySequenceLength; ++i ) {
 159             const char singleBase = Constants::BAM_DNA_LOOKUP[ ( (seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
 160             QueryBases.append(1, singleBase);
 161         }
 162     }
 163
 164     // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
 165     Qualities.clear();
 166     if ( hasQualData ) {
 167         Qualities.reserve(SupportData.QuerySequenceLength);
 168         for ( size_t i = 0; i < SupportData.QuerySequenceLength; ++i ) {
 169             const char singleQuality = static_cast<const char>(qualData[i]+33);
 170             Qualities.append(1, singleQuality);
 171         }
 172     }
 173
 174     // clear previous AlignedBases
 175     AlignedBases.clear();
 176
 177     // if QueryBases has data, build AlignedBases using CIGAR data
 178     // otherwise, AlignedBases will remain empty (this case IS allowed)
 179     if ( !QueryBases.empty() ) {
 180
 181         // resize AlignedBases
 182         AlignedBases.reserve(SupportData.QuerySequenceLength);
 183
 184         // iterate over CigarOps
 185         int k = 0;
 186         vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 187         vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 188         for ( ; cigarIter != cigarEnd; ++cigarIter ) {
 189             const CigarOp& op = (*cigarIter);
 190
 191             switch ( op.Type ) {
 192
 193                 // for 'M', 'I', '=', 'X' - write bases
 194                 case (Constants::BAM_CIGAR_MATCH_CHAR)    :
 195                 case (Constants::BAM_CIGAR_INS_CHAR)      :
 196                 case (Constants::BAM_CIGAR_SEQMATCH_CHAR) :
 197                 case (Constants::BAM_CIGAR_MISMATCH_CHAR) :
 198                     AlignedBases.append(QueryBases.substr(k, op.Length));
 199                     // fall through
 200
 201                 // for 'S' - soft clip, do not write bases
 202                 // but increment placeholder 'k'
 203                 case (Constants::BAM_CIGAR_SOFTCLIP_CHAR) :
 204                     k += op.Length;
 205                     break;
 206
 207                 // for 'D' - write gap character
 208                 case (Constants::BAM_CIGAR_DEL_CHAR) :
 209                     AlignedBases.append(op.Length, Constants::BAM_DNA_DEL);
 210                     break;
 211
 212                 // for 'P' - write padding character
 213                 case (Constants::BAM_CIGAR_PAD_CHAR) :
 214                     AlignedBases.append( op.Length, Constants::BAM_DNA_PAD );
 215                     break;
 216
 217                 // for 'N' - write N's, skip bases in original query sequence
 218                 case (Constants::BAM_CIGAR_REFSKIP_CHAR) :
 219                     AlignedBases.append( op.Length, Constants::BAM_DNA_N );
 220                     break;
 221
 222                 // for 'H' - hard clip, do nothing to AlignedBases, move to next op
 223                 case (Constants::BAM_CIGAR_HARDCLIP_CHAR) :
 224                     break;
 225
 226                 // invalid CIGAR op-code
 227                 default:
 228                     const string message = string("invalid CIGAR operation type: ") + op.Type;
 229                     SetErrorString("BamAlignment::BuildCharData", message);
 230                     return false;
 231             }
 232         }
 233     }
 234
 235     // save tag data
 236     TagData.clear();
 237     if ( hasTagData ) {
 238         if ( IsBigEndian ) {
 239             size_t i = 0;
 240             while ( i < tagDataLength ) {
 241
 242                 i += Constants::BAM_TAG_TAGSIZE;  // skip tag chars (e.g. "RG", "NM", etc.)
 243                 const char type = tagData[i];     // get tag type at position i
 244                 ++i;                              // move i past tag type
 245
 246                 switch (type) {
 247
 248                     case(Constants::BAM_TAG_TYPE_ASCII) :
 249                     case(Constants::BAM_TAG_TYPE_INT8)  :
 250                     case(Constants::BAM_TAG_TYPE_UINT8) :
 251                         // no endian swapping necessary for single-byte data
 252                         ++i;
 253                         break;
 254
 255                     case(Constants::BAM_TAG_TYPE_INT16)  :
 256                     case(Constants::BAM_TAG_TYPE_UINT16) :
 257                         BamTools::SwapEndian_16p(&tagData[i]);
 258                         i += sizeof(uint16_t);
 259                         break;
 260
 261                     case(Constants::BAM_TAG_TYPE_FLOAT)  :
 262                     case(Constants::BAM_TAG_TYPE_INT32)  :
 263                     case(Constants::BAM_TAG_TYPE_UINT32) :
 264                         BamTools::SwapEndian_32p(&tagData[i]);
 265                         i += sizeof(uint32_t);
 266                         break;
 267
 268                     case(Constants::BAM_TAG_TYPE_HEX) :
 269                     case(Constants::BAM_TAG_TYPE_STRING) :
 270                         // no endian swapping necessary for hex-string/string data
 271                         while ( tagData[i] )
 272                             ++i;
 273                         // increment one more for null terminator
 274                         ++i;
 275                         break;
 276
 277                     case(Constants::BAM_TAG_TYPE_ARRAY) :
 278
 279                     {
 280                         // read array type
 281                         const char arrayType = tagData[i];
 282                         ++i;
 283
 284                         // swap endian-ness of number of elements in place, then retrieve for loop
 285                         BamTools::SwapEndian_32p(&tagData[i]);
 286                         uint32_t numElements;
 287                         memcpy(&numElements, &tagData[i], sizeof(uint32_t));
 288                         i += sizeof(uint32_t);
 289
 290                         // swap endian-ness of array elements
 291                         for ( size_t j = 0; j < numElements; ++j ) {
 292                             switch (arrayType) {
 293                                 case (Constants::BAM_TAG_TYPE_INT8)  :
 294                                 case (Constants::BAM_TAG_TYPE_UINT8) :
 295                                     // no endian-swapping necessary
 296                                     ++i;
 297                                     break;
 298                                 case (Constants::BAM_TAG_TYPE_INT16)  :
 299                                 case (Constants::BAM_TAG_TYPE_UINT16) :
 300                                     BamTools::SwapEndian_16p(&tagData[i]);
 301                                     i += sizeof(uint16_t);
 302                                     break;
 303                                 case (Constants::BAM_TAG_TYPE_FLOAT)  :
 304                                 case (Constants::BAM_TAG_TYPE_INT32)  :
 305                                 case (Constants::BAM_TAG_TYPE_UINT32) :
 306                                     BamTools::SwapEndian_32p(&tagData[i]);
 307                                     i += sizeof(uint32_t);
 308                                     break;
 309                                 default:
 310                                     const string message = string("invalid binary array type: ") + arrayType;
 311                                     SetErrorString("BamAlignment::BuildCharData", message);
 312                                     return false;
 313                             }
 314                         }
 315
 316                         break;
 317                     }
 318
 319                     // invalid tag type-code
 320                     default :
 321                         const string message = string("invalid tag type: ") + type;
 322                         SetErrorString("BamAlignment::BuildCharData", message);
 323                         return false;
 324                 }
 325             }
 326         }
 327
 328         // store tagData in alignment
 329         TagData.resize(tagDataLength);
 330         memcpy((char*)(TagData.data()), tagData, tagDataLength);
 331     }
 332
 333     // clear core-only flag & return success
 334     SupportData.HasCoreOnly = false;
 335     return true;
 336 }
 337
 338 /*! \fn bool BamAlignment::FindTag(const std::string& tag, char*& pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) const
 339     \internal
 340
 341     Searches for requested tag in BAM tag data.
 342
 343     \param[in]     tag            requested 2-character tag name
 344     \param[in,out] pTagData       pointer to current position in BamAlignment::TagData
 345     \param[in]     tagDataLength  length of BamAlignment::TagData
 346     \param[in,out] numBytesParsed number of bytes parsed so far
 347
 348     \return \c true if found
 349
 350     \post If \a tag is found, \a pTagData will point to the byte where the tag data begins.
 351           \a numBytesParsed will correspond to the position in the full TagData string.
 352
 353 */
 354 bool BamAlignment::FindTag(const std::string& tag,
 355                            char*& pTagData,
 356                            const unsigned int& tagDataLength,
 357                            unsigned int& numBytesParsed) const
 358 {
 359
 360     while ( numBytesParsed < tagDataLength ) {
 361
 362         const char* pTagType        = pTagData;
 363         const char* pTagStorageType = pTagData + 2;
 364         pTagData       += 3;
 365         numBytesParsed += 3;
 366
 367         // check the current tag, return true on match
 368         if ( strncmp(pTagType, tag.c_str(), 2) == 0 )
 369             return true;
 370
 371         // get the storage class and find the next tag
 372         if ( *pTagStorageType == '\0' ) return false;
 373         if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
 374         if ( *pTagData == '\0' ) return false;
 375     }
 376
 377     // checked all tags, none match
 378     return false;
 379 }
 380
 381 /*! \fn int BamAlignment::GetEndPosition(bool usePadded = false, bool closedInterval = false) const
 382     \brief Calculates alignment end position, based on its starting position and CIGAR data.
 383
 384     \warning The position returned now represents a zero-based, HALF-OPEN interval.
 385     In previous versions of BamTools (0.x & 1.x) all intervals were treated
 386     as zero-based, CLOSED.
 387
 388     \param[in] usePadded      Allow inserted bases to affect the reported position. Default is
 389                               false, so that reported position stays synced with reference
 390                               coordinates.
 391     \param[in] closedInterval Setting this to true will return a 0-based end coordinate. Default is
 392                               false, so that his value represents a standard, half-open interval.
 393
 394     \return alignment end position
 395 */
 396 int BamAlignment::GetEndPosition(bool usePadded, bool closedInterval) const {
 397
 398     // initialize alignment end to starting position
 399     int alignEnd = Position;
 400
 401     // iterate over cigar operations
 402     vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 403     vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 404     for ( ; cigarIter != cigarEnd; ++cigarIter) {
 405         const CigarOp& op = (*cigarIter);
 406
 407         switch ( op.Type ) {
 408
 409             // increase end position on CIGAR chars [DMXN=]
 410             case Constants::BAM_CIGAR_DEL_CHAR      :
 411             case Constants::BAM_CIGAR_MATCH_CHAR    :
 412             case Constants::BAM_CIGAR_MISMATCH_CHAR :
 413             case Constants::BAM_CIGAR_REFSKIP_CHAR  :
 414             case Constants::BAM_CIGAR_SEQMATCH_CHAR :
 415                 alignEnd += op.Length;
 416                 break;
 417
 418             // increase end position on insertion, only if @usePadded is true
 419             case Constants::BAM_CIGAR_INS_CHAR :
 420                 if ( usePadded )
 421                     alignEnd += op.Length;
 422                 break;
 423
 424             // all other CIGAR chars do not affect end position
 425             default :
 426                 break;
 427         }
 428     }
 429
 430     // adjust for closedInterval, if requested
 431     if ( closedInterval )
 432         alignEnd -= 1;
 433
 434     // return result
 435     return alignEnd;
 436 }
 437
 438 /*! \fn std::string BamAlignment::GetErrorString(void) const
 439     \brief Returns a human-readable description of the last error that occurred
 440
 441     This method allows elimination of STDERR pollution. Developers of client code
 442     may choose how the messages are displayed to the user, if at all.
 443
 444     \return error description
 445 */
 446 std::string BamAlignment::GetErrorString(void) const {
 447     return ErrorString;
 448 }
 449
 450 /*! \fn bool BamAlignment::GetSoftClips(std::vector<int>& clipSizes, std::vector<int>& readPositions, std::vector<int>& genomePositions, bool usePadded = false) const
 451     \brief Identifies if an alignment has a soft clip. If so, identifies the
 452            sizes of the soft clips, as well as their positions in the read and reference.
 453
 454     \param[out] clipSizes       vector of the sizes of each soft clip in the alignment
 455     \param[out] readPositions   vector of the 0-based read locations of each soft clip in the alignment.
 456                                 These positions are basically indexes within the read, not genomic positions.
 457     \param[out] genomePositions vector of the 0-based genome locations of each soft clip in the alignment
 458     \param[in]  usePadded       inserted bases affect reported position. Default is false, so that
 459                                 reported position stays 'sync-ed' with reference coordinates.
 460
 461     \return \c true if any soft clips were found in the alignment
 462 */
 463 bool BamAlignment::GetSoftClips(vector<int>& clipSizes,
 464                                 vector<int>& readPositions,
 465                                 vector<int>& genomePositions,
 466                                 bool usePadded) const
 467 {
 468     // initialize positions & flags
 469     int refPosition  = Position;
 470     int readPosition = 0;
 471     bool softClipFound = false;
 472     bool firstCigarOp  = true;
 473
 474     // iterate over cigar operations
 475     vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
 476     vector<CigarOp>::const_iterator cigarEnd  = CigarData.end();
 477     for ( ; cigarIter != cigarEnd; ++cigarIter) {
 478         const CigarOp& op = (*cigarIter);
 479
 480         switch ( op.Type ) {
 481
 482             // increase both read & genome positions on CIGAR chars [DMXN=]
 483             case Constants::BAM_CIGAR_DEL_CHAR      :
 484             case Constants::BAM_CIGAR_MATCH_CHAR    :
 485             case Constants::BAM_CIGAR_MISMATCH_CHAR :
 486             case Constants::BAM_CIGAR_REFSKIP_CHAR  :
 487             case Constants::BAM_CIGAR_SEQMATCH_CHAR :
 488                 refPosition  += op.Length;
 489                 readPosition += op.Length;
 490                 break;
 491
 492             // increase read position on insertion, genome position only if @usePadded is true
 493             case Constants::BAM_CIGAR_INS_CHAR :
 494                 readPosition += op.Length;
 495                 if ( usePadded )
 496                     refPosition += op.Length;
 497                 break;
 498
 499             case Constants::BAM_CIGAR_SOFTCLIP_CHAR :
 500
 501                 softClipFound = true;
 502
 503                 //////////////////////////////////////////////////////////////////////////////
 504                 // if we are dealing with the *first* CIGAR operation
 505                 // for this alignment, we increment the read position so that
 506                 // the read and genome position of the clip are referring to the same base.
 507                 // For example, in the alignment below, the ref position would be 4, yet
 508                 //              the read position would be 0. Thus, to "sync" the two,
 509                 //              we need to increment the read position by the length of the
 510                 //              soft clip.
 511                 // Read:  ATCGTTTCGTCCCTGC
 512                 // Ref:   GGGATTTCGTCCCTGC
 513                 // Cigar: SSSSMMMMMMMMMMMM
 514                 //
 515                 // NOTE: This only needs to be done if the soft clip is the _first_ CIGAR op.
 516                 //////////////////////////////////////////////////////////////////////////////
 517                 if ( firstCigarOp )
 518                     readPosition += op.Length;
 519
 520                 // track the soft clip's size, read position, and genome position
 521                 clipSizes.push_back(op.Length);
 522                 readPositions.push_back(readPosition);
 523                 genomePositions.push_back(refPosition);
 524
 525             // any other CIGAR operations have no effect
 526             default :
 527                 break;
 528         }
 529
 530         // clear our "first pass" flag
 531         firstCigarOp = false;
 532     }
 533
 534     // return whether any soft clips found
 535     return softClipFound;
 536 }
 537
 538 /*! \fn bool BamAlignment::GetTagType(const std::string& tag, char& type) const
 539     \brief Retrieves the BAM tag type-code associated with requested tag name.
 540
 541     \param[in]  tag  2-character tag name
 542     \param[out] type retrieved (1-character) type-code
 543
 544     \return \c true if found
 545     \sa \samSpecURL for more details on reserved tag names, supported tag types, etc.
 546 */
 547 bool BamAlignment::GetTagType(const std::string& tag, char& type) const {
 548
 549     // skip if alignment is core-only
 550     if ( SupportData.HasCoreOnly ) {
 551         // TODO: set error string?
 552         return false;
 553     }
 554
 555     // skip if no tags present
 556     if ( TagData.empty() ) {
 557         // TODO: set error string?
 558         return false;
 559     }
 560
 561     // localize the tag data
 562     char* pTagData = (char*)TagData.data();
 563     const unsigned int tagDataLength = TagData.size();
 564     unsigned int numBytesParsed = 0;
 565
 566     // if tag not found, return failure
 567     if ( !FindTag(tag, pTagData, tagDataLength, numBytesParsed) ){
 568         // TODO: set error string?
 569         return false;
 570     }
 571
 572     // otherwise, retrieve & validate tag type code
 573     type = *(pTagData - 1);
 574     switch ( type ) {
 575         case (Constants::BAM_TAG_TYPE_ASCII)  :
 576         case (Constants::BAM_TAG_TYPE_INT8)   :
 577         case (Constants::BAM_TAG_TYPE_UINT8)  :
 578         case (Constants::BAM_TAG_TYPE_INT16)  :
 579         case (Constants::BAM_TAG_TYPE_UINT16) :
 580         case (Constants::BAM_TAG_TYPE_INT32)  :
 581         case (Constants::BAM_TAG_TYPE_UINT32) :
 582         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 583         case (Constants::BAM_TAG_TYPE_STRING) :
 584         case (Constants::BAM_TAG_TYPE_HEX)    :
 585         case (Constants::BAM_TAG_TYPE_ARRAY)  :
 586             return true;
 587
 588         // unknown tag type
 589         default:
 590             const string message = string("invalid tag type: ") + type;
 591             SetErrorString("BamAlignment::GetTagType", message);
 592             return false;
 593     }
 594 }
 595
 596 /*! \fn bool BamAlignment::HasTag(const std::string& tag) const
 597     \brief Returns true if alignment has a record for requested tag.
 598
 599     \param[in] tag 2-character tag name
 600     \return \c true if alignment has a record for tag
 601 */
 602 bool BamAlignment::HasTag(const std::string& tag) const {
 603
 604     // return false if no tag data present
 605     if ( SupportData.HasCoreOnly || TagData.empty() )
 606         return false;
 607
 608     // localize the tag data for lookup
 609     char* pTagData = (char*)TagData.data();
 610     const unsigned int tagDataLength = TagData.size();
 611     unsigned int numBytesParsed = 0;
 612
 613     // if result of tag lookup
 614     return FindTag(tag, pTagData, tagDataLength, numBytesParsed);
 615 }
 616
 617 /*! \fn bool BamAlignment::IsDuplicate(void) const
 618     \return \c true if this read is a PCR duplicate
 619 */
 620 bool BamAlignment::IsDuplicate(void) const {
 621     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_DUPLICATE) != 0 );
 622 }
 623
 624 /*! \fn bool BamAlignment::IsFailedQC(void) const
 625     \return \c true if this read failed quality control
 626 */
 627 bool BamAlignment::IsFailedQC(void) const {
 628     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_QC_FAILED) != 0 );
 629 }
 630
 631 /*! \fn bool BamAlignment::IsFirstMate(void) const
 632     \return \c true if alignment is first mate on paired-end read
 633 */
 634 bool BamAlignment::IsFirstMate(void) const {
 635     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_1) != 0 );
 636 }
 637
 638 /*! \fn bool BamAlignment::IsMapped(void) const
 639     \return \c true if alignment is mapped
 640 */
 641 bool BamAlignment::IsMapped(void) const {
 642     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_UNMAPPED) == 0 );
 643 }
 644
 645 /*! \fn bool BamAlignment::IsMateMapped(void) const
 646     \return \c true if alignment's mate is mapped
 647 */
 648 bool BamAlignment::IsMateMapped(void) const {
 649     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_UNMAPPED) == 0 );
 650 }
 651
 652 /*! \fn bool BamAlignment::IsMateReverseStrand(void) const
 653     \return \c true if alignment's mate mapped to reverse strand
 654 */
 655 bool BamAlignment::IsMateReverseStrand(void) const {
 656     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND) != 0 );
 657 }
 658
 659 /*! \fn bool BamAlignment::IsPaired(void) const
 660     \return \c true if alignment part of paired-end read
 661 */
 662 bool BamAlignment::IsPaired(void) const {
 663     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PAIRED) != 0 );
 664 }
 665
 666 /*! \fn bool BamAlignment::IsPrimaryAlignment(void) const
 667     \return \c true if reported position is primary alignment
 668 */
 669 bool BamAlignment::IsPrimaryAlignment(void) const  {
 670     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_SECONDARY) == 0 );
 671 }
 672
 673 /*! \fn bool BamAlignment::IsProperPair(void) const
 674     \return \c true if alignment is part of read that satisfied paired-end resolution
 675 */
 676 bool BamAlignment::IsProperPair(void) const {
 677     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PROPER_PAIR) != 0 );
 678 }
 679
 680 /*! \fn bool BamAlignment::IsReverseStrand(void) const
 681     \return \c true if alignment mapped to reverse strand
 682 */
 683 bool BamAlignment::IsReverseStrand(void) const {
 684     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_REVERSE_STRAND) != 0 );
 685 }
 686
 687 /*! \fn bool BamAlignment::IsSecondMate(void) const
 688     \return \c true if alignment is second mate on read
 689 */
 690 bool BamAlignment::IsSecondMate(void) const {
 691     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_2) != 0 );
 692 }
 693
 694 /*! \fn bool BamAlignment::IsValidSize(const std::string& tag, const std::string& type) const
 695     \internal
 696
 697     Checks that tag name & type strings are expected sizes.
 698
 699     \param tag[in]  BAM tag name
 700     \param type[in] BAM tag type-code
 701     \return \c true if both input strings are valid sizes
 702 */
 703 bool BamAlignment::IsValidSize(const std::string& tag, const std::string& type) const {
 704     return (tag.size()  == Constants::BAM_TAG_TAGSIZE) &&
 705            (type.size() == Constants::BAM_TAG_TYPESIZE);
 706 }
 707
 708 /*! \fn void BamAlignment::RemoveTag(const std::string& tag)
 709     \brief Removes field from BAM tags.
 710
 711     \param[in] tag 2-character name of field to remove
 712 */
 713 void BamAlignment::RemoveTag(const std::string& tag) {
 714
 715     // if char data not populated, do that first
 716     if ( SupportData.HasCoreOnly )
 717         BuildCharData();
 718
 719     // skip if no tags available
 720     if ( TagData.empty() )
 721         return;
 722
 723     // localize the tag data
 724     char* pOriginalTagData = (char*)TagData.data();
 725     char* pTagData = pOriginalTagData;
 726     const unsigned int originalTagDataLength = TagData.size();
 727     unsigned int newTagDataLength = 0;
 728     unsigned int numBytesParsed = 0;
 729
 730     // skip if tag not found
 731     if  ( !FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) )
 732         return;
 733
 734     // otherwise, remove it
 735     RaiiBuffer newTagData(originalTagDataLength);
 736
 737     // copy original tag data up til desired tag
 738     pTagData       -= 3;
 739     numBytesParsed -= 3;
 740     const unsigned int beginningTagDataLength = numBytesParsed;
 741     newTagDataLength += beginningTagDataLength;
 742     memcpy(newTagData.Buffer, pOriginalTagData, numBytesParsed);
 743
 744     // attemp to skip to next tag
 745     const char* pTagStorageType = pTagData + 2;
 746     pTagData       += 3;
 747     numBytesParsed += 3;
 748     if ( SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) {
 749
 750         // squeeze remaining tag data
 751         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 752         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 753         memcpy(newTagData.Buffer + beginningTagDataLength, pTagData, endTagDataLength );
 754
 755         // save modified tag data in alignment
 756         TagData.assign(newTagData.Buffer, beginningTagDataLength + endTagDataLength);
 757     }
 758 }
 759
 760 /*! \fn void BamAlignment::SetErrorString(const std::string& where, const std::string& what) const
 761     \internal
 762
 763     Sets a formatted error string for this alignment.
 764
 765     \param[in] where class/method where error occurred
 766     \param[in] what  description of error
 767 */
 768 void BamAlignment::SetErrorString(const std::string& where, const std::string& what) const {
 769     static const string SEPARATOR = ": ";
 770     ErrorString = where + SEPARATOR + what;
 771 }
 772
 773 /*! \fn void BamAlignment::SetIsDuplicate(bool ok)
 774     \brief Sets value of "PCR duplicate" flag to \a ok.
 775 */
 776 void BamAlignment::SetIsDuplicate(bool ok) {
 777     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_DUPLICATE;
 778     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_DUPLICATE;
 779 }
 780
 781 /*! \fn void BamAlignment::SetIsFailedQC(bool ok)
 782     \brief Sets "failed quality control" flag to \a ok.
 783 */
 784 void BamAlignment::SetIsFailedQC(bool ok) {
 785     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_QC_FAILED;
 786     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_QC_FAILED;
 787 }
 788
 789 /*! \fn void BamAlignment::SetIsFirstMate(bool ok)
 790     \brief Sets "alignment is first mate" flag to \a ok.
 791 */
 792 void BamAlignment::SetIsFirstMate(bool ok) {
 793     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_1;
 794     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_1;
 795 }
 796
 797 /*! \fn void BamAlignment::SetIsMapped(bool ok)
 798     \brief Sets "alignment is mapped" flag to \a ok.
 799 */
 800 void BamAlignment::SetIsMapped(bool ok) {
 801     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_UNMAPPED;
 802     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_UNMAPPED;
 803 }
 804
 805 /*! \fn void BamAlignment::SetIsMateMapped(bool ok)
 806     \brief Sets "alignment's mate is mapped" flag to \a ok.
 807 */
 808 void BamAlignment::SetIsMateMapped(bool ok) {
 809     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
 810     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
 811 }
 812
 813 /*! \fn void BamAlignment::SetIsMateReverseStrand(bool ok)
 814     \brief Sets "alignment's mate mapped to reverse strand" flag to \a ok.
 815 */
 816 void BamAlignment::SetIsMateReverseStrand(bool ok) {
 817     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
 818     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
 819 }
 820
 821 /*! \fn void BamAlignment::SetIsPaired(bool ok)
 822     \brief Sets "alignment part of paired-end read" flag to \a ok.
 823 */
 824 void BamAlignment::SetIsPaired(bool ok) {
 825     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PAIRED;
 826     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PAIRED;
 827 }
 828
 829 /*! \fn void BamAlignment::SetIsPrimaryAlignment(bool ok)
 830     \brief Sets "position is primary alignment" flag to \a ok.
 831 */
 832 void BamAlignment::SetIsPrimaryAlignment(bool ok) {
 833     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_SECONDARY;
 834     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_SECONDARY;
 835 }
 836
 837 /*! \fn void BamAlignment::SetIsProperPair(bool ok)
 838     \brief Sets "alignment is part of read that satisfied paired-end resolution" flag to \a ok.
 839 */
 840 void BamAlignment::SetIsProperPair(bool ok) {
 841     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PROPER_PAIR;
 842     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PROPER_PAIR;
 843 }
 844
 845 /*! \fn void BamAlignment::SetIsReverseStrand(bool ok)
 846     \brief Sets "alignment mapped to reverse strand" flag to \a ok.
 847 */
 848 void BamAlignment::SetIsReverseStrand(bool ok) {
 849     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_REVERSE_STRAND;
 850     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_REVERSE_STRAND;
 851 }
 852
 853 /*! \fn void BamAlignment::SetIsSecondMate(bool ok)
 854     \brief Sets "alignment is second mate on read" flag to \a ok.
 855 */
 856 void BamAlignment::SetIsSecondMate(bool ok) {
 857     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_2;
 858     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_2;
 859 }
 860
 861 /*! \fn bool BamAlignment::SkipToNextTag(const char storageType, char*& pTagData, unsigned int& numBytesParsed) const
 862     \internal
 863
 864     Moves to next available tag in tag data string
 865
 866     \param[in]     storageType    BAM tag type-code that determines how far to move cursor
 867     \param[in,out] pTagData       pointer to current position (cursor) in tag string
 868     \param[in,out] numBytesParsed report of how many bytes were parsed (cumulatively)
 869
 870     \return \c if storageType was a recognized BAM tag type
 871
 872     \post \a pTagData       will point to the byte where the next tag data begins.
 873           \a numBytesParsed will correspond to the cursor's position in the full TagData string.
 874 */
 875 bool BamAlignment::SkipToNextTag(const char storageType,
 876                                  char*& pTagData,
 877                                  unsigned int& numBytesParsed) const
 878 {
 879     switch (storageType) {
 880
 881         case (Constants::BAM_TAG_TYPE_ASCII) :
 882         case (Constants::BAM_TAG_TYPE_INT8)  :
 883         case (Constants::BAM_TAG_TYPE_UINT8) :
 884             ++numBytesParsed;
 885             ++pTagData;
 886             break;
 887
 888         case (Constants::BAM_TAG_TYPE_INT16)  :
 889         case (Constants::BAM_TAG_TYPE_UINT16) :
 890             numBytesParsed += sizeof(uint16_t);
 891             pTagData       += sizeof(uint16_t);
 892             break;
 893
 894         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 895         case (Constants::BAM_TAG_TYPE_INT32)  :
 896         case (Constants::BAM_TAG_TYPE_UINT32) :
 897             numBytesParsed += sizeof(uint32_t);
 898             pTagData       += sizeof(uint32_t);
 899             break;
 900
 901         case (Constants::BAM_TAG_TYPE_STRING) :
 902         case (Constants::BAM_TAG_TYPE_HEX)    :
 903             while( *pTagData ) {
 904                 ++numBytesParsed;
 905                 ++pTagData;
 906             }
 907             // increment for null-terminator
 908             ++numBytesParsed;
 909             ++pTagData;
 910             break;
 911
 912         case (Constants::BAM_TAG_TYPE_ARRAY) :
 913
 914         {
 915             // read array type
 916             const char arrayType = *pTagData;
 917             ++numBytesParsed;
 918             ++pTagData;
 919
 920             // read number of elements
 921             int32_t numElements;
 922             memcpy(&numElements, pTagData, sizeof(uint32_t)); // already endian-swapped, if needed
 923             numBytesParsed += sizeof(uint32_t);
 924             pTagData       += sizeof(uint32_t);
 925
 926             // calculate number of bytes to skip
 927             int bytesToSkip = 0;
 928             switch (arrayType) {
 929                 case (Constants::BAM_TAG_TYPE_INT8)  :
 930                 case (Constants::BAM_TAG_TYPE_UINT8) :
 931                     bytesToSkip = numElements;
 932                     break;
 933                 case (Constants::BAM_TAG_TYPE_INT16)  :
 934                 case (Constants::BAM_TAG_TYPE_UINT16) :
 935                     bytesToSkip = numElements*sizeof(uint16_t);
 936                     break;
 937                 case (Constants::BAM_TAG_TYPE_FLOAT)  :
 938                 case (Constants::BAM_TAG_TYPE_INT32)  :
 939                 case (Constants::BAM_TAG_TYPE_UINT32) :
 940                     bytesToSkip = numElements*sizeof(uint32_t);
 941                     break;
 942                 default:
 943                     const string message = string("invalid binary array type: ") + arrayType;
 944                     SetErrorString("BamAlignment::SkipToNextTag", message);
 945                     return false;
 946             }
 947
 948             // skip binary array contents
 949             numBytesParsed += bytesToSkip;
 950             pTagData       += bytesToSkip;
 951             break;
 952         }
 953
 954         default:
 955             const string message = string("invalid tag type: ") + storageType;
 956             SetErrorString("BamAlignment::SkipToNextTag", message);
 957             return false;
 958     }
 959
 960     // if we get here, tag skipped OK - return success
 961     return true;
 962 }