src/api/BamAlignment.h

   1 // ***************************************************************************
   2 // BamAlignment.h (c) 2009 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 4 October 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Provides the BamAlignment data structure
   8 // ***************************************************************************
   9
  10 #ifndef BAMALIGNMENT_H
  11 #define BAMALIGNMENT_H
  12
  13 #include <api/api_global.h>
  14 #include <api/BamAux.h>
  15 #include <api/BamConstants.h>
  16 #include <cctype>
  17 #include <cstdio>
  18 #include <cstdlib>
  19 #include <cstring>
  20 #include <iostream>
  21 #include <string>
  22 #include <vector>
  23
  24 namespace BamTools {
  25
  26 // forward declaration of BamAlignment's friend classes
  27 namespace Internal {
  28     class BamReaderPrivate;
  29     class BamWriterPrivate;
  30 } // namespace Internal
  31
  32 // BamAlignment data structure
  33 struct API_EXPORT BamAlignment {
  34
  35     // constructors & destructor
  36     public:
  37         BamAlignment(void);
  38         BamAlignment(const BamAlignment& other);
  39         ~BamAlignment(void);
  40
  41     // queries against alignment flags
  42     public:
  43         bool IsDuplicate(void) const;           // returns true if this read is a PCR duplicate
  44         bool IsFailedQC(void) const;            // returns true if this read failed quality control
  45         bool IsFirstMate(void) const;           // returns true if alignment is first mate on read
  46         bool IsMapped(void) const;              // returns true if alignment is mapped
  47         bool IsMateMapped(void) const;          // returns true if alignment's mate is mapped
  48         bool IsMateReverseStrand(void) const;   // returns true if alignment's mate mapped to reverse strand
  49         bool IsPaired(void) const;              // returns true if alignment part of paired-end read
  50         bool IsPrimaryAlignment(void) const;    // returns true if reported position is primary alignment
  51         bool IsProperPair(void) const;          // returns true if alignment is part of read that satisfied paired-end resolution
  52         bool IsReverseStrand(void) const;       // returns true if alignment mapped to reverse strand
  53         bool IsSecondMate(void) const;          // returns true if alignment is second mate on read
  54
  55     // manipulate alignment flags
  56     public:
  57         void SetIsDuplicate(bool ok);           // sets value of "PCR duplicate" flag
  58         void SetIsFailedQC(bool ok);            // sets value of "failed quality control" flag
  59         void SetIsFirstMate(bool ok);           // sets value of "alignment is first mate" flag
  60         void SetIsMapped(bool ok);              // sets value of "alignment is mapped" flag
  61         void SetIsMateMapped(bool ok);          // sets value of "alignment's mate is mapped" flag
  62         void SetIsMateReverseStrand(bool ok);   // sets value of "alignment's mate mapped to reverse strand" flag
  63         void SetIsPaired(bool ok);              // sets value of "alignment part of paired-end read" flag
  64         void SetIsPrimaryAlignment(bool ok);    // sets value of "position is primary alignment" flag
  65         void SetIsProperPair(bool ok);          // sets value of "alignment is part of read that satisfied paired-end resolution" flag
  66         void SetIsReverseStrand(bool ok);       // sets value of "alignment mapped to reverse strand" flag
  67         void SetIsSecondMate(bool ok);          // sets value of "alignment is second mate on read" flag
  68
  69         // legacy methods (consider deprecated, but still available)
  70         void SetIsMateUnmapped(bool ok);        // complement of using SetIsMateMapped()
  71         void SetIsSecondaryAlignment(bool ok);  // complement of using SetIsPrimaryAlignment()
  72         void SetIsUnmapped(bool ok);            // complement of using SetIsMapped()
  73
  74     // tag data access methods
  75     public:
  76
  77         // add a new tag
  78         template<typename T> bool AddTag(const std::string& tag,
  79                                          const std::string& type,
  80                                          const T& value);
  81         template<typename T> bool AddTag(const std::string& tag,
  82                                          const std::vector<T>& values);
  83
  84         // edit (or append) tag
  85         template<typename T> bool EditTag(const std::string& tag,
  86                                           const std::string& type,
  87                                           const T& value);
  88         template<typename T> bool EditTag(const std::string& tag,
  89                                           const std::vector<T>& values);
  90
  91         // retrieves tag data
  92         template<typename T> bool GetTag(const std::string& tag,
  93                                          T& destination) const;
  94         template<typename T> bool GetTag(const std::string& tag,
  95                                          std::vector<T>& destination) const;
  96
  97         // retrieves the BAM type-code for requested tag
  98         // (returns whether or not tag exists, and type-code is valid)
  99         bool GetTagType(const std::string& tag, char& type) const;
 100
 101         // legacy methods (consider deprecated, but still available)
 102         bool GetEditDistance(uint32_t& editDistance) const;         // retrieves value of "NM" tag
 103         bool GetReadGroup(std::string& readGroup) const;            // retrieves value of "RG" tag
 104
 105         // returns true if alignment has a record for this tag name
 106         bool HasTag(const std::string& tag) const;
 107
 108         // removes a tag
 109         bool RemoveTag(const std::string& tag);
 110
 111     // additional methods
 112     public:
 113         // populates alignment string fields
 114         bool BuildCharData(void);
 115
 116         // calculates alignment end position
 117         int GetEndPosition(bool usePadded = false, bool zeroBased = true) const;
 118
 119     // public data fields
 120     public:
 121         std::string Name;               // read name
 122         int32_t     Length;             // length of query sequence
 123         std::string QueryBases;         // 'original' sequence (as reported from sequencing machine)
 124         std::string AlignedBases;       // 'aligned' sequence (includes any indels, padding, clipping)
 125         std::string Qualities;          // FASTQ qualities (ASCII characters, not numeric values)
 126         std::string TagData;            // tag data (use provided methods to query/modify)
 127         int32_t     RefID;              // ID number for reference sequence
 128         int32_t     Position;           // position (0-based) where alignment starts
 129         uint16_t    Bin;                // BAM (standard) index bin number for this alignment
 130         uint16_t    MapQuality;         // mapping quality score
 131         uint32_t    AlignmentFlag;      // alignment bit-flag (use provided methods to query/modify)
 132         std::vector<CigarOp> CigarData; // CIGAR operations for this alignment
 133         int32_t     MateRefID;          // ID number for reference sequence where alignment's mate was aligned
 134         int32_t     MatePosition;       // position (0-based) where alignment's mate starts
 135         int32_t     InsertSize;         // mate-pair insert size
 136         std::string Filename;           // name of BAM file which this alignment comes from
 137
 138     //! \cond
 139     // internal utility methods
 140     private:
 141         static bool FindTag(const std::string& tag,
 142                             char*& pTagData,
 143                             const unsigned int& tagDataLength,
 144                             unsigned int& numBytesParsed);
 145         static bool IsValidSize(const std::string& tag,
 146                                 const std::string& type);
 147         static bool SkipToNextTag(const char storageType,
 148                            char*& pTagData,
 149                            unsigned int& numBytesParsed);
 150
 151     // internal data
 152     private:
 153
 154         struct BamAlignmentSupportData {
 155
 156             // data members
 157             std::string AllCharData;
 158             uint32_t    BlockLength;
 159             uint32_t    NumCigarOperations;
 160             uint32_t    QueryNameLength;
 161             uint32_t    QuerySequenceLength;
 162             bool        HasCoreOnly;
 163
 164             // constructor
 165             BamAlignmentSupportData(void)
 166                 : BlockLength(0)
 167                 , NumCigarOperations(0)
 168                 , QueryNameLength(0)
 169                 , QuerySequenceLength(0)
 170                 , HasCoreOnly(false)
 171             { }
 172         };
 173         BamAlignmentSupportData SupportData;
 174         friend class Internal::BamReaderPrivate;
 175         friend class Internal::BamWriterPrivate;
 176     //! \endcond
 177 };
 178
 179 // ---------------------------------------------------------
 180 // BamAlignment tag access methods
 181
 182 template<typename T>
 183 inline bool BamAlignment::AddTag(const std::string& tag,
 184                                  const std::string& type,
 185                                  const T& value)
 186 {
 187     // if char data not populated, do that first
 188     if ( SupportData.HasCoreOnly )
 189         BuildCharData();
 190
 191     // validate tag/type size & that storage type code is OK for T
 192     if ( !IsValidSize(tag, type) ) return false;
 193     if ( !TagTypeHelper<T>::CanConvertTo(type.at(0)) )
 194         return false;
 195
 196     // localize the tag data
 197     char* pTagData = (char*)TagData.data();
 198     const unsigned int tagDataLength = TagData.size();
 199     unsigned int numBytesParsed = 0;
 200
 201     // if tag already exists, return false
 202     // use EditTag explicitly instead
 203     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 204         return false;
 205
 206     // otherwise, convert value to string
 207     union { T value; char valueBuffer[sizeof(T)]; } un;
 208     un.value = value;
 209
 210     // copy original tag data to temp buffer
 211     const std::string newTag = tag + type;
 212     const int newTagDataLength = tagDataLength + newTag.size() + sizeof(T); // leave room for new T
 213     char* originalTagData = new char[newTagDataLength];
 214     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 215
 216     // append newTag
 217     strcat(originalTagData + tagDataLength, newTag.data());
 218     memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(T));
 219
 220     // store temp buffer back in TagData
 221     const char* newTagData = (const char*)originalTagData;
 222     TagData.assign(newTagData, newTagDataLength);
 223
 224     // clean up & return success
 225     delete[] originalTagData;
 226     return true;
 227 }
 228
 229 template<>
 230 inline bool BamAlignment::AddTag<std::string>(const std::string& tag,
 231                                               const std::string& type,
 232                                               const std::string& value)
 233 {
 234     // if char data not populated, do that first
 235     if ( SupportData.HasCoreOnly )
 236         BuildCharData();
 237
 238     // validate tag/type size & that storage type code is OK for string
 239     if ( !IsValidSize(tag, type) ) return false;
 240     if ( !TagTypeHelper<std::string>::CanConvertTo(type.at(0)) )
 241         return false;
 242
 243     // localize the tag data
 244     char* pTagData = (char*)TagData.data();
 245     const unsigned int tagDataLength = TagData.size();
 246     unsigned int numBytesParsed = 0;
 247
 248     // if tag already exists, return false
 249     // use EditTag explicitly instead
 250     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 251         return false;
 252
 253     // otherwise, copy tag data to temp buffer
 254     const std::string newTag = tag + type + value;
 255     const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
 256     char* originalTagData = new char[newTagDataLength];
 257     memcpy(originalTagData, TagData.c_str(), tagDataLength + 1);    // '+1' for TagData null-term
 258
 259     // append newTag
 260     strcat(originalTagData + tagDataLength, newTag.data());  // removes original null-term, appends newTag + null-term
 261
 262     // store temp buffer back in TagData
 263     const char* newTagData = (const char*)originalTagData;
 264     TagData.assign(newTagData, newTagDataLength);
 265
 266     // clean up & return success
 267     delete[] originalTagData;
 268     return true;
 269 }
 270
 271 template<typename T>
 272 inline bool BamAlignment::AddTag(const std::string& tag,
 273                                  const std::vector<T>& values)
 274 {
 275     // if char data not populated, do that first
 276     if ( SupportData.HasCoreOnly )
 277         BuildCharData();
 278
 279     // check for valid tag name length
 280     if ( tag.size() != Constants::BAM_TAG_TAGSIZE )
 281         return false;
 282
 283     // localize the tag data
 284     char* pTagData = (char*)TagData.data();
 285     const unsigned int tagDataLength = TagData.size();
 286     unsigned int numBytesParsed = 0;
 287
 288     // if tag already exists, return false
 289     // use EditTag explicitly instead
 290     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 291         return false;
 292
 293     // build new tag's base information
 294     char newTagBase[Constants::BAM_TAG_ARRAYBASE_SIZE];
 295     memcpy( newTagBase, tag.c_str(), Constants::BAM_TAG_TAGSIZE );
 296     newTagBase[2] = Constants::BAM_TAG_TYPE_ARRAY;
 297     newTagBase[3] = TagTypeHelper<T>::TypeCode();
 298
 299     // add number of array elements to newTagBase
 300     const int32_t numElements  = values.size();
 301     memcpy(newTagBase + 4, &numElements, sizeof(int32_t));
 302
 303     // copy current TagData string to temp buffer, leaving room for new tag's contents
 304     const int newTagDataLength = tagDataLength +
 305                                  Constants::BAM_TAG_ARRAYBASE_SIZE +
 306                                  numElements*sizeof(T);
 307     char* originalTagData = new char[newTagDataLength];
 308     memcpy(originalTagData, TagData.c_str(), tagDataLength+1); // '+1' for TagData's null-term
 309
 310     // write newTagBase (removes old null term)
 311     strcat(originalTagData + tagDataLength, (const char*)newTagBase);
 312
 313     // add vector elements to tag
 314     int elementsBeginOffset = tagDataLength + Constants::BAM_TAG_ARRAYBASE_SIZE;
 315     for ( int i = 0 ; i < numElements; ++i ) {
 316         const T& value = values.at(i);
 317         memcpy(originalTagData + elementsBeginOffset + i*sizeof(T), &value, sizeof(T));
 318     }
 319
 320     // store temp buffer back in TagData
 321     const char* newTagData = (const char*)originalTagData;
 322     TagData.assign(newTagData, newTagDataLength);
 323
 324     // cleanup & return success
 325     delete[] originalTagData;
 326     return true;
 327 }
 328
 329 template<typename T>
 330 inline bool BamAlignment::EditTag(const std::string& tag,
 331                                   const std::string& type,
 332                                   const T& value)
 333 {
 334     // if char data not populated, do that first
 335     if ( SupportData.HasCoreOnly )
 336         BuildCharData();
 337
 338     // remove existing tag if present, then append tag with new value
 339     if ( HasTag(tag) )
 340         RemoveTag(tag);
 341     return AddTag(tag, type, value);
 342 }
 343
 344 template<typename T>
 345 inline bool BamAlignment::EditTag(const std::string& tag,
 346                                   const std::vector<T>& values)
 347 {
 348     // if char data not populated, do that first
 349     if ( SupportData.HasCoreOnly )
 350         BuildCharData();
 351
 352     // remove existing tag if present, then append tag with new values
 353     if ( HasTag(tag) )
 354         RemoveTag(tag);
 355     return AddTag(tag, values);
 356 }
 357
 358 template<typename T>
 359 inline bool BamAlignment::GetTag(const std::string& tag,
 360                                  T& destination) const
 361 {
 362     // skip if core-only or no tags present
 363     if ( SupportData.HasCoreOnly || TagData.empty() )
 364         return false;
 365
 366     // localize the tag data
 367     char* pTagData = (char*)TagData.data();
 368     const unsigned int tagDataLength = TagData.size();
 369     unsigned int numBytesParsed = 0;
 370
 371     // return failure if tag not found
 372     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 373         return false;
 374
 375     // otherwise try to copy data into destination
 376     const char type = *(pTagData - 1);
 377     if ( !TagTypeHelper<T>::CanConvertFrom(type) )
 378         return false;
 379     int destinationLength = 0;
 380     switch ( type ) {
 381
 382         // 1 byte data
 383         case (Constants::BAM_TAG_TYPE_ASCII) :
 384         case (Constants::BAM_TAG_TYPE_INT8)  :
 385         case (Constants::BAM_TAG_TYPE_UINT8) :
 386             destinationLength = 1;
 387             break;
 388
 389         // 2 byte data
 390         case (Constants::BAM_TAG_TYPE_INT16)  :
 391         case (Constants::BAM_TAG_TYPE_UINT16) :
 392             destinationLength = 2;
 393             break;
 394
 395         // 4 byte data
 396         case (Constants::BAM_TAG_TYPE_INT32)  :
 397         case (Constants::BAM_TAG_TYPE_UINT32) :
 398         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 399             destinationLength = 4;
 400             break;
 401
 402         // var-length types not supported for numeric destination
 403         case (Constants::BAM_TAG_TYPE_STRING) :
 404         case (Constants::BAM_TAG_TYPE_HEX)    :
 405         case (Constants::BAM_TAG_TYPE_ARRAY)  :
 406             std::cerr << "BamAlignment ERROR: cannot store tag of type " << type
 407                  << " in integer destination" << std::endl;
 408             return false;
 409
 410         // unrecognized tag type
 411         default:
 412             std::cerr << "BamAlignment ERROR: unknown tag type encountered: "
 413                  << type << std::endl;
 414             return false;
 415     }
 416
 417     // store in destination
 418     destination = 0;
 419     memcpy(&destination, pTagData, destinationLength);
 420
 421     // return success
 422     return true;
 423 }
 424
 425 template<>
 426 inline bool BamAlignment::GetTag<std::string>(const std::string& tag,
 427                                               std::string& destination) const
 428 {
 429     // skip if core-only or no tags present
 430     if ( SupportData.HasCoreOnly || TagData.empty() )
 431         return false;
 432
 433     // localize the tag data
 434     char* pTagData = (char*)TagData.data();
 435     const unsigned int tagDataLength = TagData.size();
 436     unsigned int numBytesParsed = 0;
 437
 438     // return failure if tag not found
 439     if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 440         return false;
 441
 442     // otherwise copy data into destination
 443     const unsigned int dataLength = strlen(pTagData);
 444     destination.clear();
 445     destination.resize(dataLength);
 446     memcpy( (char*)destination.data(), pTagData, dataLength );
 447
 448     // return success
 449     return true;
 450 }
 451
 452 // retrieves "binary-array" tag data
 453 template<typename T>
 454 inline bool BamAlignment::GetTag(const std::string& tag,
 455                                  std::vector<T>& destination) const
 456 {
 457     // skip if core-only or no tags present
 458     if ( SupportData.HasCoreOnly || TagData.empty() )
 459         return false;
 460
 461     // localize the tag data
 462     char* pTagData = (char*)TagData.data();
 463     const unsigned int tagDataLength = TagData.size();
 464     unsigned int numBytesParsed = 0;
 465
 466     // return false if tag not found
 467     if ( !FindTag(tag, pTagData, tagDataLength, numBytesParsed) )
 468         return false;
 469
 470     // check that tag is array type
 471     const char tagType = *(pTagData - 1);
 472     if ( tagType != Constants::BAM_TAG_TYPE_ARRAY ) {
 473         std::cerr << "BamAlignment ERROR: Cannot store non-array data from tag: "
 474              << tag << " in array destination" << std::endl;
 475         return false;
 476     }
 477
 478     // calculate length of each element in tag's array
 479     const char elementType = *pTagData;
 480     if ( !TagTypeHelper<T>::CanConvertFrom(elementType) )
 481         return false;
 482     ++pTagData;
 483     int elementLength = 0;
 484     switch ( elementType ) {
 485         case (Constants::BAM_TAG_TYPE_ASCII) :
 486         case (Constants::BAM_TAG_TYPE_INT8)  :
 487         case (Constants::BAM_TAG_TYPE_UINT8) :
 488             elementLength = sizeof(uint8_t);
 489             break;
 490
 491         case (Constants::BAM_TAG_TYPE_INT16)  :
 492         case (Constants::BAM_TAG_TYPE_UINT16) :
 493             elementLength = sizeof(uint16_t);
 494             break;
 495
 496         case (Constants::BAM_TAG_TYPE_INT32)  :
 497         case (Constants::BAM_TAG_TYPE_UINT32) :
 498         case (Constants::BAM_TAG_TYPE_FLOAT)  :
 499             elementLength = sizeof(uint32_t);
 500             break;
 501
 502         // var-length types not supported for numeric destination
 503         case (Constants::BAM_TAG_TYPE_STRING) :
 504         case (Constants::BAM_TAG_TYPE_HEX)    :
 505         case (Constants::BAM_TAG_TYPE_ARRAY)  :
 506             std::cerr << "BamAlignment ERROR: array element type: " << elementType
 507                  << " cannot be stored in integer value" << std::endl;
 508             return false;
 509
 510         // unknown tag type
 511         default:
 512             std::cerr << "BamAlignment ERROR: unknown element type encountered: "
 513                  << elementType << std::endl;
 514             return false;
 515     }
 516
 517     // get number of elements
 518     int32_t numElements;
 519     memcpy(&numElements, pTagData, sizeof(int32_t));
 520     pTagData += 4;
 521     destination.clear();
 522     destination.reserve(numElements);
 523
 524     // read in elements
 525     T value;
 526     for ( int i = 0 ; i < numElements; ++i ) {
 527         memcpy(&value, pTagData, sizeof(T));
 528         pTagData += sizeof(T);
 529         destination.push_back(value);
 530     }
 531
 532     // return success
 533     return true;
 534 }
 535
 536 typedef std::vector<BamAlignment> BamAlignmentVector;
 537
 538 } // namespace BamTools
 539
 540 #endif // BAMALIGNMENT_H