src/api/BamAlignment.cpp

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