src/api/BamAlignment.cpp

   1 // ***************************************************************************
   2 // BamAlignment.cpp (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 4 April 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 bool BamAlignment::GetTagType(const std::string& tag, char& type) const
 555     \brief Retrieves the BAM tag type-code associated with requested tag name.
 556
 557     \param[in]  tag  2-character tag name
 558     \param[out] type retrieved (1-character) type-code
 559
 560     \return \c true if found
 561     \sa \samSpecURL for more details on reserved tag names, supported tag types, etc.
 562 */
 563 bool BamAlignment::GetTagType(const std::string& tag, char& type) const {
 564
 565     // skip if alignment is core-only
 566     if ( SupportData.HasCoreOnly ) {
 567         // TODO: set error string?
 568         return false;
 569     }
 570
 571     // skip if no tags present
 572     if ( TagData.empty() ) {
 573         // TODO: set error string?
 574         return false;
 575     }
 576
 577     // localize the tag data
 578     char* pTagData = (char*)TagData.data();
 579     const unsigned int tagDataLength = TagData.size();
 580     unsigned int numBytesParsed = 0;
 581
 582     // if tag not found, return failure
 583     if ( !FindTag(tag, pTagData, tagDataLength, numBytesParsed) ){
 584         // TODO: set error string?
 585         return false;
 586     }
 587
 588     // otherwise, retrieve & validate tag type code
 589     type = *(pTagData - 1);
 590     switch ( type ) {
 591         case (Constants::BAM_TAG_TYPE_ASCII)  :
 592         case (Constants::BAM_TAG_TYPE_INT8)   :
 593         case (Constants::BAM_TAG_TYPE_UINT8)  :
 594         case (Constants::BAM_TAG_TYPE_INT16)  :
 595         case (Constants::BAM_TAG_TYPE_UINT16) :
 596         case (Constants::BAM_TAG_TYPE_INT32)  :
 597         case (Constants::BAM_TAG_TYPE_UINT32) :
 598         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 599         case (Constants::BAM_TAG_TYPE_STRING) :
 600         case (Constants::BAM_TAG_TYPE_HEX)    :
 601         case (Constants::BAM_TAG_TYPE_ARRAY)  :
 602             return true;
 603
 604         // unknown tag type
 605         default:
 606             const string message = string("invalid tag type: ") + type;
 607             SetErrorString("BamAlignment::GetTagType", message);
 608             return false;
 609     }
 610 }
 611
 612 /*! \fn bool BamAlignment::HasTag(const std::string& tag) const
 613     \brief Returns true if alignment has a record for requested tag.
 614
 615     \param[in] tag 2-character tag name
 616     \return \c true if alignment has a record for tag
 617 */
 618 bool BamAlignment::HasTag(const std::string& tag) const {
 619
 620     // return false if no tag data present
 621     if ( SupportData.HasCoreOnly || TagData.empty() )
 622         return false;
 623
 624     // localize the tag data for lookup
 625     char* pTagData = (char*)TagData.data();
 626     const unsigned int tagDataLength = TagData.size();
 627     unsigned int numBytesParsed = 0;
 628
 629     // if result of tag lookup
 630     return FindTag(tag, pTagData, tagDataLength, numBytesParsed);
 631 }
 632
 633 /*! \fn bool BamAlignment::IsDuplicate(void) const
 634     \return \c true if this read is a PCR duplicate
 635 */
 636 bool BamAlignment::IsDuplicate(void) const {
 637     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_DUPLICATE) != 0 );
 638 }
 639
 640 /*! \fn bool BamAlignment::IsFailedQC(void) const
 641     \return \c true if this read failed quality control
 642 */
 643 bool BamAlignment::IsFailedQC(void) const {
 644     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_QC_FAILED) != 0 );
 645 }
 646
 647 /*! \fn bool BamAlignment::IsFirstMate(void) const
 648     \return \c true if alignment is first mate on paired-end read
 649 */
 650 bool BamAlignment::IsFirstMate(void) const {
 651     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_1) != 0 );
 652 }
 653
 654 /*! \fn bool BamAlignment::IsMapped(void) const
 655     \return \c true if alignment is mapped
 656 */
 657 bool BamAlignment::IsMapped(void) const {
 658     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_UNMAPPED) == 0 );
 659 }
 660
 661 /*! \fn bool BamAlignment::IsMateMapped(void) const
 662     \return \c true if alignment's mate is mapped
 663 */
 664 bool BamAlignment::IsMateMapped(void) const {
 665     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_UNMAPPED) == 0 );
 666 }
 667
 668 /*! \fn bool BamAlignment::IsMateReverseStrand(void) const
 669     \return \c true if alignment's mate mapped to reverse strand
 670 */
 671 bool BamAlignment::IsMateReverseStrand(void) const {
 672     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND) != 0 );
 673 }
 674
 675 /*! \fn bool BamAlignment::IsPaired(void) const
 676     \return \c true if alignment part of paired-end read
 677 */
 678 bool BamAlignment::IsPaired(void) const {
 679     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PAIRED) != 0 );
 680 }
 681
 682 /*! \fn bool BamAlignment::IsPrimaryAlignment(void) const
 683     \return \c true if reported position is primary alignment
 684 */
 685 bool BamAlignment::IsPrimaryAlignment(void) const  {
 686     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_SECONDARY) == 0 );
 687 }
 688
 689 /*! \fn bool BamAlignment::IsProperPair(void) const
 690     \return \c true if alignment is part of read that satisfied paired-end resolution
 691 */
 692 bool BamAlignment::IsProperPair(void) const {
 693     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_PROPER_PAIR) != 0 );
 694 }
 695
 696 /*! \fn bool BamAlignment::IsReverseStrand(void) const
 697     \return \c true if alignment mapped to reverse strand
 698 */
 699 bool BamAlignment::IsReverseStrand(void) const {
 700     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_REVERSE_STRAND) != 0 );
 701 }
 702
 703 /*! \fn bool BamAlignment::IsSecondMate(void) const
 704     \return \c true if alignment is second mate on read
 705 */
 706 bool BamAlignment::IsSecondMate(void) const {
 707     return ( (AlignmentFlag & Constants::BAM_ALIGNMENT_READ_2) != 0 );
 708 }
 709
 710 /*! \fn bool BamAlignment::IsValidSize(const std::string& tag, const std::string& type) const
 711     \internal
 712
 713     Checks that tag name & type strings are expected sizes.
 714
 715     \param tag[in]  BAM tag name
 716     \param type[in] BAM tag type-code
 717     \return \c true if both input strings are valid sizes
 718 */
 719 bool BamAlignment::IsValidSize(const std::string& tag, const std::string& type) const {
 720     return (tag.size()  == Constants::BAM_TAG_TAGSIZE) &&
 721            (type.size() == Constants::BAM_TAG_TYPESIZE);
 722 }
 723
 724 /*! \fn void BamAlignment::RemoveTag(const std::string& tag)
 725     \brief Removes field from BAM tags.
 726
 727     \param[in] tag 2-character name of field to remove
 728 */
 729 void BamAlignment::RemoveTag(const std::string& tag) {
 730
 731     // if char data not populated, do that first
 732     if ( SupportData.HasCoreOnly )
 733         BuildCharData();
 734
 735     // skip if no tags available
 736     if ( TagData.empty() )
 737         return;
 738
 739     // localize the tag data
 740     char* pOriginalTagData = (char*)TagData.data();
 741     char* pTagData = pOriginalTagData;
 742     const unsigned int originalTagDataLength = TagData.size();
 743     unsigned int newTagDataLength = 0;
 744     unsigned int numBytesParsed = 0;
 745
 746     // skip if tag not found
 747     if  ( !FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) )
 748         return;
 749
 750     // otherwise, remove it
 751     RaiiBuffer newTagData(originalTagDataLength);
 752
 753     // copy original tag data up til desired tag
 754     pTagData       -= 3;
 755     numBytesParsed -= 3;
 756     const unsigned int beginningTagDataLength = numBytesParsed;
 757     newTagDataLength += beginningTagDataLength;
 758     memcpy(newTagData.Buffer, pOriginalTagData, numBytesParsed);
 759
 760     // attemp to skip to next tag
 761     const char* pTagStorageType = pTagData + 2;
 762     pTagData       += 3;
 763     numBytesParsed += 3;
 764     if ( SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) {
 765
 766         // squeeze remaining tag data
 767         const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
 768         const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
 769         memcpy(newTagData.Buffer + beginningTagDataLength, pTagData, endTagDataLength );
 770
 771         // save modified tag data in alignment
 772         TagData.assign(newTagData.Buffer, beginningTagDataLength + endTagDataLength);
 773     }
 774 }
 775
 776 /*! \fn void BamAlignment::SetErrorString(const std::string& where, const std::string& what) const
 777     \internal
 778
 779     Sets a formatted error string for this alignment.
 780
 781     \param[in] where class/method where error occurred
 782     \param[in] what  description of error
 783 */
 784 void BamAlignment::SetErrorString(const std::string& where, const std::string& what) const {
 785     static const string SEPARATOR = ": ";
 786     ErrorString = where + SEPARATOR + what;
 787 }
 788
 789 /*! \fn void BamAlignment::SetIsDuplicate(bool ok)
 790     \brief Sets value of "PCR duplicate" flag to \a ok.
 791 */
 792 void BamAlignment::SetIsDuplicate(bool ok) {
 793     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_DUPLICATE;
 794     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_DUPLICATE;
 795 }
 796
 797 /*! \fn void BamAlignment::SetIsFailedQC(bool ok)
 798     \brief Sets "failed quality control" flag to \a ok.
 799 */
 800 void BamAlignment::SetIsFailedQC(bool ok) {
 801     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_QC_FAILED;
 802     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_QC_FAILED;
 803 }
 804
 805 /*! \fn void BamAlignment::SetIsFirstMate(bool ok)
 806     \brief Sets "alignment is first mate" flag to \a ok.
 807 */
 808 void BamAlignment::SetIsFirstMate(bool ok) {
 809     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_1;
 810     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_1;
 811 }
 812
 813 /*! \fn void BamAlignment::SetIsMapped(bool ok)
 814     \brief Sets "alignment is mapped" flag to \a ok.
 815 */
 816 void BamAlignment::SetIsMapped(bool ok) {
 817     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_UNMAPPED;
 818     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_UNMAPPED;
 819 }
 820
 821 /*! \fn void BamAlignment::SetIsMateMapped(bool ok)
 822     \brief Sets "alignment's mate is mapped" flag to \a ok.
 823 */
 824 void BamAlignment::SetIsMateMapped(bool ok) {
 825     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
 826     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_UNMAPPED;
 827 }
 828
 829 /*! \fn void BamAlignment::SetIsMateReverseStrand(bool ok)
 830     \brief Sets "alignment's mate mapped to reverse strand" flag to \a ok.
 831 */
 832 void BamAlignment::SetIsMateReverseStrand(bool ok) {
 833     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
 834     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_MATE_REVERSE_STRAND;
 835 }
 836
 837 /*! \fn void BamAlignment::SetIsPaired(bool ok)
 838     \brief Sets "alignment part of paired-end read" flag to \a ok.
 839 */
 840 void BamAlignment::SetIsPaired(bool ok) {
 841     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PAIRED;
 842     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PAIRED;
 843 }
 844
 845 /*! \fn void BamAlignment::SetIsPrimaryAlignment(bool ok)
 846     \brief Sets "position is primary alignment" flag to \a ok.
 847 */
 848 void BamAlignment::SetIsPrimaryAlignment(bool ok) {
 849     if (ok) AlignmentFlag &= ~Constants::BAM_ALIGNMENT_SECONDARY;
 850     else    AlignmentFlag |=  Constants::BAM_ALIGNMENT_SECONDARY;
 851 }
 852
 853 /*! \fn void BamAlignment::SetIsProperPair(bool ok)
 854     \brief Sets "alignment is part of read that satisfied paired-end resolution" flag to \a ok.
 855 */
 856 void BamAlignment::SetIsProperPair(bool ok) {
 857     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_PROPER_PAIR;
 858     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_PROPER_PAIR;
 859 }
 860
 861 /*! \fn void BamAlignment::SetIsReverseStrand(bool ok)
 862     \brief Sets "alignment mapped to reverse strand" flag to \a ok.
 863 */
 864 void BamAlignment::SetIsReverseStrand(bool ok) {
 865     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_REVERSE_STRAND;
 866     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_REVERSE_STRAND;
 867 }
 868
 869 /*! \fn void BamAlignment::SetIsSecondMate(bool ok)
 870     \brief Sets "alignment is second mate on read" flag to \a ok.
 871 */
 872 void BamAlignment::SetIsSecondMate(bool ok) {
 873     if (ok) AlignmentFlag |=  Constants::BAM_ALIGNMENT_READ_2;
 874     else    AlignmentFlag &= ~Constants::BAM_ALIGNMENT_READ_2;
 875 }
 876
 877 /*! \fn bool BamAlignment::SkipToNextTag(const char storageType, char*& pTagData, unsigned int& numBytesParsed) const
 878     \internal
 879
 880     Moves to next available tag in tag data string
 881
 882     \param[in]     storageType    BAM tag type-code that determines how far to move cursor
 883     \param[in,out] pTagData       pointer to current position (cursor) in tag string
 884     \param[in,out] numBytesParsed report of how many bytes were parsed (cumulatively)
 885
 886     \return \c if storageType was a recognized BAM tag type
 887
 888     \post \a pTagData       will point to the byte where the next tag data begins.
 889           \a numBytesParsed will correspond to the cursor's position in the full TagData string.
 890 */
 891 bool BamAlignment::SkipToNextTag(const char storageType,
 892                                  char*& pTagData,
 893                                  unsigned int& numBytesParsed) const
 894 {
 895     switch (storageType) {
 896
 897         case (Constants::BAM_TAG_TYPE_ASCII) :
 898         case (Constants::BAM_TAG_TYPE_INT8)  :
 899         case (Constants::BAM_TAG_TYPE_UINT8) :
 900             ++numBytesParsed;
 901             ++pTagData;
 902             break;
 903
 904         case (Constants::BAM_TAG_TYPE_INT16)  :
 905         case (Constants::BAM_TAG_TYPE_UINT16) :
 906             numBytesParsed += sizeof(uint16_t);
 907             pTagData       += sizeof(uint16_t);
 908             break;
 909
 910         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 911         case (Constants::BAM_TAG_TYPE_INT32)  :
 912         case (Constants::BAM_TAG_TYPE_UINT32) :
 913             numBytesParsed += sizeof(uint32_t);
 914             pTagData       += sizeof(uint32_t);
 915             break;
 916
 917         case (Constants::BAM_TAG_TYPE_STRING) :
 918         case (Constants::BAM_TAG_TYPE_HEX)    :
 919             while( *pTagData ) {
 920                 ++numBytesParsed;
 921                 ++pTagData;
 922             }
 923             // increment for null-terminator
 924             ++numBytesParsed;
 925             ++pTagData;
 926             break;
 927
 928         case (Constants::BAM_TAG_TYPE_ARRAY) :
 929
 930         {
 931             // read array type
 932             const char arrayType = *pTagData;
 933             ++numBytesParsed;
 934             ++pTagData;
 935
 936             // read number of elements
 937             int32_t numElements;
 938             memcpy(&numElements, pTagData, sizeof(uint32_t)); // already endian-swapped, if needed
 939             numBytesParsed += sizeof(uint32_t);
 940             pTagData       += sizeof(uint32_t);
 941
 942             // calculate number of bytes to skip
 943             int bytesToSkip = 0;
 944             switch (arrayType) {
 945                 case (Constants::BAM_TAG_TYPE_INT8)  :
 946                 case (Constants::BAM_TAG_TYPE_UINT8) :
 947                     bytesToSkip = numElements;
 948                     break;
 949                 case (Constants::BAM_TAG_TYPE_INT16)  :
 950                 case (Constants::BAM_TAG_TYPE_UINT16) :
 951                     bytesToSkip = numElements*sizeof(uint16_t);
 952                     break;
 953                 case (Constants::BAM_TAG_TYPE_FLOAT)  :
 954                 case (Constants::BAM_TAG_TYPE_INT32)  :
 955                 case (Constants::BAM_TAG_TYPE_UINT32) :
 956                     bytesToSkip = numElements*sizeof(uint32_t);
 957                     break;
 958                 default:
 959                     const string message = string("invalid binary array type: ") + arrayType;
 960                     SetErrorString("BamAlignment::SkipToNextTag", message);
 961                     return false;
 962             }
 963
 964             // skip binary array contents
 965             numBytesParsed += bytesToSkip;
 966             pTagData       += bytesToSkip;
 967             break;
 968         }
 969
 970         default:
 971             const string message = string("invalid tag type: ") + storageType;
 972             SetErrorString("BamAlignment::SkipToNextTag", message);
 973             return false;
 974     }
 975
 976     // if we get here, tag skipped OK - return success
 977     return true;
 978 }