src/api/BamAlignment.cpp

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