1 // ***************************************************************************
\r
2 // BamAux.h (c) 2009 Derek Barnett, Michael Str�mberg
\r
3 // Marth Lab, Department of Biology, Boston College
\r
4 // All rights reserved.
\r
5 // ---------------------------------------------------------------------------
\r
6 // Last modified: 16 September 2010 (DB)
\r
7 // ---------------------------------------------------------------------------
\r
8 // Provides the basic constants, data structures, etc. for using BAM files
\r
9 // ***************************************************************************
\r
21 #include <exception>
\r
29 // Platform-specific type definitions
\r
30 #ifndef BAMTOOLS_TYPES
\r
31 #define BAMTOOLS_TYPES
\r
33 typedef char int8_t;
\r
34 typedef unsigned char uint8_t;
\r
35 typedef short int16_t;
\r
36 typedef unsigned short uint16_t;
\r
37 typedef int int32_t;
\r
38 typedef unsigned int uint32_t;
\r
39 typedef long long int64_t;
\r
40 typedef unsigned long long uint64_t;
\r
44 #endif // BAMTOOLS_TYPES
\r
46 namespace BamTools {
\r
49 const int BAM_CIGAR_SHIFT = 4;
\r
50 const int BAM_CIGAR_MASK = ((1 << BAM_CIGAR_SHIFT) - 1);
\r
51 const int BAM_CORE_SIZE = 32;
\r
55 struct BamAlignment {
\r
57 // constructors & destructor
\r
60 BamAlignment(const BamAlignment& other);
\r
61 ~BamAlignment(void);
\r
63 // Queries against alignment flags
\r
65 bool IsDuplicate(void) const; // Returns true if this read is a PCR duplicate
\r
66 bool IsFailedQC(void) const; // Returns true if this read failed quality control
\r
67 bool IsFirstMate(void) const; // Returns true if alignment is first mate on read
\r
68 bool IsMapped(void) const; // Returns true if alignment is mapped
\r
69 bool IsMateMapped(void) const; // Returns true if alignment's mate is mapped
\r
70 bool IsMateReverseStrand(void) const; // Returns true if alignment's mate mapped to reverse strand
\r
71 bool IsPaired(void) const; // Returns true if alignment part of paired-end read
\r
72 bool IsPrimaryAlignment(void) const; // Returns true if reported position is primary alignment
\r
73 bool IsProperPair(void) const; // Returns true if alignment is part of read that satisfied paired-end resolution
\r
74 bool IsReverseStrand(void) const; // Returns true if alignment mapped to reverse strand
\r
75 bool IsSecondMate(void) const; // Returns true if alignment is second mate on read
\r
77 // Manipulate alignment flags
\r
79 void SetIsDuplicate(bool ok); // Sets "PCR duplicate" flag
\r
80 void SetIsFailedQC(bool ok); // Sets "failed quality control" flag
\r
81 void SetIsFirstMate(bool ok); // Sets "alignment is first mate" flag
\r
82 void SetIsMateUnmapped(bool ok); // Sets "alignment's mate is mapped" flag
\r
83 void SetIsMateReverseStrand(bool ok); // Sets "alignment's mate mapped to reverse strand" flag
\r
84 void SetIsPaired(bool ok); // Sets "alignment part of paired-end read" flag
\r
85 void SetIsProperPair(bool ok); // Sets "alignment is part of read that satisfied paired-end resolution" flag
\r
86 void SetIsReverseStrand(bool ok); // Sets "alignment mapped to reverse strand" flag
\r
87 void SetIsSecondaryAlignment(bool ok); // Sets "position is primary alignment" flag
\r
88 void SetIsSecondMate(bool ok); // Sets "alignment is second mate on read" flag
\r
89 void SetIsUnmapped(bool ok); // Sets "alignment is mapped" flag
\r
91 // Tag data access methods
\r
93 // -------------------------------------------------------------------------------------
\r
94 // N.B. - The following tag-modifying methods may not be used on BamAlignments fetched
\r
95 // using BamReader::GetNextAlignmentCore(). Attempting to use them will not result in
\r
96 // error message (to keep output clean) but will ALWAYS return false. Only user-
\r
97 // generated BamAlignments or those retrieved using BamReader::GetNextAlignment() are valid.
\r
99 // add tag data (create new TAG entry with TYPE and VALUE)
\r
100 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
101 // returns true if new data added, false if error or TAG already exists
\r
102 // N.B. - will NOT modify existing tag. Use EditTag() instead
\r
103 bool AddTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
104 bool AddTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
105 bool AddTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
106 bool AddTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
108 // edit tag data (sets existing TAG with TYPE to VALUE or adds new TAG if not already present)
\r
109 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
110 // returns true if edit was successfaul, false if error
\r
111 bool EditTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
112 bool EditTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
113 bool EditTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
114 bool EditTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
116 // specific tag data access methods - these only remain for legacy support
\r
117 bool GetEditDistance(uint32_t& editDistance) const; // get "NM" tag data (implemented as GetTag("NM", editDistance))
\r
118 bool GetReadGroup(std::string& readGroup) const; // get "RG" tag data (implemented as GetTag("RG", readGroup))
\r
120 // generic tag data access methods
\r
121 bool GetTag(const std::string& tag, std::string& destination) const; // access variable-length char or hex strings
\r
122 bool GetTag(const std::string& tag, uint32_t& destination) const; // access unsigned integer data
\r
123 bool GetTag(const std::string& tag, int32_t& destination) const; // access signed integer data
\r
124 bool GetTag(const std::string& tag, float& destination) const; // access floating point data
\r
127 // returns true if removal was successful, false if error
\r
128 // N.B. - returns false if TAG does not exist (no removal can occur)
\r
129 bool RemoveTag(const std::string& tag);
\r
131 // Additional data access methods
\r
133 // calculates alignment end position, based on starting position and CIGAR operations
\r
134 // @zeroBased - if true, returns 0-based coordinate; else returns 1-based
\r
135 int GetEndPosition(bool usePadded = false, bool zeroBased = true) const;
\r
137 // 'internal' utility methods
\r
139 static bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed);
\r
140 static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);
\r
144 std::string Name; // Read name
\r
145 int32_t Length; // Query length
\r
146 std::string QueryBases; // 'Original' sequence (as reported from sequencing machine)
\r
147 std::string AlignedBases; // 'Aligned' sequence (includes any indels, padding, clipping)
\r
148 std::string Qualities; // FASTQ qualities (ASCII characters, not numeric values)
\r
149 std::string TagData; // Tag data (accessor methods will pull the requested information out)
\r
150 int32_t RefID; // ID number for reference sequence
\r
151 int32_t Position; // Position (0-based) where alignment starts
\r
152 uint16_t Bin; // Bin in BAM file where this alignment resides
\r
153 uint16_t MapQuality; // Mapping quality score
\r
154 uint32_t AlignmentFlag; // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate
\r
155 std::vector<CigarOp> CigarData; // CIGAR operations for this alignment
\r
156 int32_t MateRefID; // ID number for reference sequence where alignment's mate was aligned
\r
157 int32_t MatePosition; // Position (0-based) where alignment's mate starts
\r
158 int32_t InsertSize; // Mate-pair insert size
\r
162 struct BamAlignmentSupportData {
\r
165 std::string AllCharData;
\r
166 uint32_t BlockLength;
\r
167 uint32_t NumCigarOperations;
\r
168 uint32_t QueryNameLength;
\r
169 uint32_t QuerySequenceLength;
\r
173 BamAlignmentSupportData(void)
\r
175 , NumCigarOperations(0)
\r
176 , QueryNameLength(0)
\r
177 , QuerySequenceLength(0)
\r
178 , HasCoreOnly(false)
\r
182 // contains raw character data & lengths
\r
183 BamAlignmentSupportData SupportData;
\r
185 // allow these classes access to BamAlignment private members (SupportData)
\r
186 // but client code should not need to touch this data
\r
187 friend class BamReader;
\r
188 friend class BamWriter;
\r
190 // Alignment flag query constants
\r
191 // Use the get/set methods above instead
\r
196 , MATE_UNMAPPED = 8
\r
198 , MATE_REVERSE = 32
\r
203 , DUPLICATE = 1024
\r
207 // ----------------------------------------------------------------
\r
208 // Auxiliary data structs & typedefs
\r
213 char Type; // Operation type (MIDNSHP)
\r
214 uint32_t Length; // Operation length (number of bases)
\r
217 CigarOp(const char type = '\0',
\r
218 const uint32_t length = 0)
\r
227 std::string RefName; // Name of reference sequence
\r
228 int32_t RefLength; // Length of reference sequence
\r
229 bool RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence
\r
232 RefData(const int32_t& length = 0,
\r
234 : RefLength(length)
\r
235 , RefHasAlignments(ok)
\r
239 typedef std::vector<RefData> RefVector;
\r
240 typedef std::vector<BamAlignment> BamAlignmentVector;
\r
251 BamRegion(const int& leftID = -1,
\r
252 const int& leftPos = -1,
\r
253 const int& rightID = -1,
\r
254 const int& rightPos = -1)
\r
255 : LeftRefID(leftID)
\r
256 , LeftPosition(leftPos)
\r
257 , RightRefID(rightID)
\r
258 , RightPosition(rightPos)
\r
261 // member functions
\r
262 void clear(void) { LeftRefID = -1; LeftPosition = -1; RightRefID = -1; RightPosition = -1; }
\r
263 bool isLeftBoundSpecified(void) const { return ( LeftRefID != -1 && LeftPosition != -1 ); }
\r
264 bool isNull(void) const { return ( !isLeftBoundSpecified() && !isRightBoundSpecified() ); }
\r
265 bool isRightBoundSpecified(void) const { return ( RightRefID != -1 && RightPosition != -1 ); }
\r
268 // ----------------------------------------------------------------
\r
269 // Added: 3-35-2010 DWB
\r
270 // Fixed: Routines to provide endian-correctness
\r
271 // ----------------------------------------------------------------
\r
273 // returns true if system is big endian
\r
274 inline bool SystemIsBigEndian(void) {
\r
275 const uint16_t one = 0x0001;
\r
276 return ((*(char*) &one) == 0 );
\r
279 // swaps endianness of 16-bit value 'in place'
\r
280 inline void SwapEndian_16(int16_t& x) {
\r
281 x = ((x >> 8) | (x << 8));
\r
284 inline void SwapEndian_16(uint16_t& x) {
\r
285 x = ((x >> 8) | (x << 8));
\r
288 // swaps endianness of 32-bit value 'in-place'
\r
289 inline void SwapEndian_32(int32_t& x) {
\r
291 ((x << 8) & 0x00FF0000) |
\r
292 ((x >> 8) & 0x0000FF00) |
\r
297 inline void SwapEndian_32(uint32_t& x) {
\r
299 ((x << 8) & 0x00FF0000) |
\r
300 ((x >> 8) & 0x0000FF00) |
\r
305 // swaps endianness of 64-bit value 'in-place'
\r
306 inline void SwapEndian_64(int64_t& x) {
\r
308 ((x << 40) & 0x00FF000000000000ll) |
\r
309 ((x << 24) & 0x0000FF0000000000ll) |
\r
310 ((x << 8) & 0x000000FF00000000ll) |
\r
311 ((x >> 8) & 0x00000000FF000000ll) |
\r
312 ((x >> 24) & 0x0000000000FF0000ll) |
\r
313 ((x >> 40) & 0x000000000000FF00ll) |
\r
318 inline void SwapEndian_64(uint64_t& x) {
\r
320 ((x << 40) & 0x00FF000000000000ll) |
\r
321 ((x << 24) & 0x0000FF0000000000ll) |
\r
322 ((x << 8) & 0x000000FF00000000ll) |
\r
323 ((x >> 8) & 0x00000000FF000000ll) |
\r
324 ((x >> 24) & 0x0000000000FF0000ll) |
\r
325 ((x >> 40) & 0x000000000000FF00ll) |
\r
330 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)
\r
331 inline void SwapEndian_16p(char* data) {
\r
332 uint16_t& value = (uint16_t&)*data;
\r
333 SwapEndian_16(value);
\r
336 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)
\r
337 inline void SwapEndian_32p(char* data) {
\r
338 uint32_t& value = (uint32_t&)*data;
\r
339 SwapEndian_32(value);
\r
342 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)
\r
343 inline void SwapEndian_64p(char* data) {
\r
344 uint64_t& value = (uint64_t&)*data;
\r
345 SwapEndian_64(value);
\r
348 inline bool FileExists(const std::string& filename) {
\r
349 std::ifstream f(filename.c_str(), std::ifstream::in);
\r
353 // ----------------------------------------------------------------
\r
354 // BamAlignment member methods
\r
356 // constructors & destructor
\r
357 inline BamAlignment::BamAlignment(void) { }
\r
359 inline BamAlignment::BamAlignment(const BamAlignment& other)
\r
361 , Length(other.Length)
\r
362 , QueryBases(other.QueryBases)
\r
363 , AlignedBases(other.AlignedBases)
\r
364 , Qualities(other.Qualities)
\r
365 , TagData(other.TagData)
\r
366 , RefID(other.RefID)
\r
367 , Position(other.Position)
\r
369 , MapQuality(other.MapQuality)
\r
370 , AlignmentFlag(other.AlignmentFlag)
\r
371 , CigarData(other.CigarData)
\r
372 , MateRefID(other.MateRefID)
\r
373 , MatePosition(other.MatePosition)
\r
374 , InsertSize(other.InsertSize)
\r
375 , SupportData(other.SupportData)
\r
378 inline BamAlignment::~BamAlignment(void) { }
\r
380 // Queries against alignment flags
\r
381 inline bool BamAlignment::IsDuplicate(void) const { return ( (AlignmentFlag & DUPLICATE) != 0 ); }
\r
382 inline bool BamAlignment::IsFailedQC(void) const { return ( (AlignmentFlag & QC_FAILED) != 0 ); }
\r
383 inline bool BamAlignment::IsFirstMate(void) const { return ( (AlignmentFlag & READ_1) != 0 ); }
\r
384 inline bool BamAlignment::IsMapped(void) const { return ( (AlignmentFlag & UNMAPPED) == 0 ); }
\r
385 inline bool BamAlignment::IsMateMapped(void) const { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }
\r
386 inline bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE) != 0 ); }
\r
387 inline bool BamAlignment::IsPaired(void) const { return ( (AlignmentFlag & PAIRED) != 0 ); }
\r
388 inline bool BamAlignment::IsPrimaryAlignment(void) const { return ( (AlignmentFlag & SECONDARY) == 0 ); }
\r
389 inline bool BamAlignment::IsProperPair(void) const { return ( (AlignmentFlag & PROPER_PAIR) != 0 ); }
\r
390 inline bool BamAlignment::IsReverseStrand(void) const { return ( (AlignmentFlag & REVERSE) != 0 ); }
\r
391 inline bool BamAlignment::IsSecondMate(void) const { return ( (AlignmentFlag & READ_2) != 0 ); }
\r
393 // Manipulate alignment flags
\r
394 inline void BamAlignment::SetIsDuplicate(bool ok) { if (ok) AlignmentFlag |= DUPLICATE; else AlignmentFlag &= ~DUPLICATE; }
\r
395 inline void BamAlignment::SetIsFailedQC(bool ok) { if (ok) AlignmentFlag |= QC_FAILED; else AlignmentFlag &= ~QC_FAILED; }
\r
396 inline void BamAlignment::SetIsFirstMate(bool ok) { if (ok) AlignmentFlag |= READ_1; else AlignmentFlag &= ~READ_1; }
\r
397 inline void BamAlignment::SetIsMateUnmapped(bool ok) { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }
\r
398 inline void BamAlignment::SetIsMateReverseStrand(bool ok) { if (ok) AlignmentFlag |= MATE_REVERSE; else AlignmentFlag &= ~MATE_REVERSE; }
\r
399 inline void BamAlignment::SetIsPaired(bool ok) { if (ok) AlignmentFlag |= PAIRED; else AlignmentFlag &= ~PAIRED; }
\r
400 inline void BamAlignment::SetIsProperPair(bool ok) { if (ok) AlignmentFlag |= PROPER_PAIR; else AlignmentFlag &= ~PROPER_PAIR; }
\r
401 inline void BamAlignment::SetIsReverseStrand(bool ok) { if (ok) AlignmentFlag |= REVERSE; else AlignmentFlag &= ~REVERSE; }
\r
402 inline void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY; else AlignmentFlag &= ~SECONDARY; }
\r
403 inline void BamAlignment::SetIsSecondMate(bool ok) { if (ok) AlignmentFlag |= READ_2; else AlignmentFlag &= ~READ_2; }
\r
404 inline void BamAlignment::SetIsUnmapped(bool ok) { if (ok) AlignmentFlag |= UNMAPPED; else AlignmentFlag &= ~UNMAPPED; }
\r
406 // calculates alignment end position, based on starting position and CIGAR operations
\r
408 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {
\r
410 // initialize alignment end to starting position
\r
411 int alignEnd = Position;
\r
413 // iterate over cigar operations
\r
414 std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
\r
415 std::vector<CigarOp>::const_iterator cigarEnd = CigarData.end();
\r
416 for ( ; cigarIter != cigarEnd; ++cigarIter) {
\r
417 const char cigarType = (*cigarIter).Type;
\r
418 if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {
\r
419 alignEnd += (*cigarIter).Length;
\r
421 else if ( usePadded && cigarType == 'I' ) {
\r
422 alignEnd += (*cigarIter).Length;
\r
426 // adjust for zeroBased, if necessary
\r
428 return alignEnd - 1;
\r
434 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
436 if ( SupportData.HasCoreOnly ) return false;
\r
437 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
438 if ( type != "Z" && type != "H" ) return false;
\r
440 // localize the tag data
\r
441 char* pTagData = (char*)TagData.data();
\r
442 const unsigned int tagDataLength = TagData.size();
\r
443 unsigned int numBytesParsed = 0;
\r
445 // if tag already exists, return false
\r
446 // use EditTag explicitly instead
\r
447 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
449 // otherwise, copy tag data to temp buffer
\r
450 std::string newTag = tag + type + value;
\r
451 const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
\r
452 char originalTagData[newTagDataLength];
\r
453 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
456 strcat(originalTagData + tagDataLength, newTag.data()); // removes original null-term, appends newTag + null-term
\r
458 // store temp buffer back in TagData
\r
459 const char* newTagData = (const char*)originalTagData;
\r
460 TagData.assign(newTagData, newTagDataLength);
\r
467 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
469 if ( SupportData.HasCoreOnly ) return false;
\r
470 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
471 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
473 // localize the tag data
\r
474 char* pTagData = (char*)TagData.data();
\r
475 const unsigned int tagDataLength = TagData.size();
\r
476 unsigned int numBytesParsed = 0;
\r
478 // if tag already exists, return false
\r
479 // use EditTag explicitly instead
\r
480 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
482 // otherwise, convert value to string
\r
483 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
486 // copy original tag data to temp buffer
\r
487 std::string newTag = tag + type;
\r
488 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
\r
489 char originalTagData[newTagDataLength];
\r
490 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
493 strcat(originalTagData + tagDataLength, newTag.data());
\r
494 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
\r
496 // store temp buffer back in TagData
\r
497 const char* newTagData = (const char*)originalTagData;
\r
498 TagData.assign(newTagData, newTagDataLength);
\r
505 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
506 return AddTag(tag, type, (const uint32_t&)value);
\r
510 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
\r
512 if ( SupportData.HasCoreOnly ) return false;
\r
513 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
514 if ( type == "Z" || type == "H" ) return false;
\r
516 // localize the tag data
\r
517 char* pTagData = (char*)TagData.data();
\r
518 const unsigned int tagDataLength = TagData.size();
\r
519 unsigned int numBytesParsed = 0;
\r
521 // if tag already exists, return false
\r
522 // use EditTag explicitly instead
\r
523 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
525 // otherwise, convert value to string
\r
526 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
529 // copy original tag data to temp buffer
\r
530 std::string newTag = tag + type;
\r
531 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
\r
532 char originalTagData[newTagDataLength];
\r
533 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
536 strcat(originalTagData + tagDataLength, newTag.data());
\r
537 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
\r
539 // store temp buffer back in TagData
\r
540 const char* newTagData = (const char*)originalTagData;
\r
541 TagData.assign(newTagData, newTagDataLength);
\r
548 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
550 if ( SupportData.HasCoreOnly ) return false;
\r
551 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
552 if ( type != "Z" && type != "H" ) return false;
\r
554 // localize the tag data
\r
555 char* pOriginalTagData = (char*)TagData.data();
\r
556 char* pTagData = pOriginalTagData;
\r
557 const unsigned int originalTagDataLength = TagData.size();
\r
559 unsigned int newTagDataLength = 0;
\r
560 unsigned int numBytesParsed = 0;
\r
562 // if tag found, store data in readGroup, return success
\r
563 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
565 // make sure array is more than big enough
\r
566 char newTagData[originalTagDataLength + value.size()];
\r
568 // copy original tag data up til desired tag
\r
569 const unsigned int beginningTagDataLength = numBytesParsed;
\r
570 newTagDataLength += beginningTagDataLength;
\r
571 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
573 // copy new VALUE in place of current tag data
\r
574 const unsigned int dataLength = strlen(value.c_str());
\r
575 memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
\r
577 // skip to next tag (if tag for removal is last, return true)
\r
578 const char* pTagStorageType = pTagData - 1;
\r
579 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
581 // copy everything from current tag (the next one after tag for removal) to end
\r
582 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
583 const unsigned int endTagOffset = beginningTagDataLength + dataLength + 1;
\r
584 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
585 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
587 // ensure null-terminator
\r
588 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
590 // save new tag data
\r
591 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
595 // tag not found, attempt AddTag
\r
596 else return AddTag(tag, type, value);
\r
600 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
602 if ( SupportData.HasCoreOnly ) return false;
\r
603 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
604 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
606 // localize the tag data
\r
607 char* pOriginalTagData = (char*)TagData.data();
\r
608 char* pTagData = pOriginalTagData;
\r
609 const unsigned int originalTagDataLength = TagData.size();
\r
611 unsigned int newTagDataLength = 0;
\r
612 unsigned int numBytesParsed = 0;
\r
614 // if tag found, store data in readGroup, return success
\r
615 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
617 // make sure array is more than big enough
\r
618 char newTagData[originalTagDataLength + sizeof(value)];
\r
620 // copy original tag data up til desired tag
\r
621 const unsigned int beginningTagDataLength = numBytesParsed;
\r
622 newTagDataLength += beginningTagDataLength;
\r
623 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
625 // copy new VALUE in place of current tag data
\r
626 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
628 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
\r
630 // skip to next tag (if tag for removal is last, return true)
\r
631 const char* pTagStorageType = pTagData - 1;
\r
632 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
634 // copy everything from current tag (the next one after tag for removal) to end
\r
635 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
636 const unsigned int endTagOffset = beginningTagDataLength + sizeof(unsigned int);
\r
637 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
638 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
640 // ensure null-terminator
\r
641 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
643 // save new tag data
\r
644 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
648 // tag not found, attempt AddTag
\r
649 else return AddTag(tag, type, value);
\r
653 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
654 return EditTag(tag, type, (const uint32_t&)value);
\r
658 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
\r
660 if ( SupportData.HasCoreOnly ) return false;
\r
661 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
662 if ( type == "Z" || type == "H" ) return false;
\r
664 // localize the tag data
\r
665 char* pOriginalTagData = (char*)TagData.data();
\r
666 char* pTagData = pOriginalTagData;
\r
667 const unsigned int originalTagDataLength = TagData.size();
\r
669 unsigned int newTagDataLength = 0;
\r
670 unsigned int numBytesParsed = 0;
\r
672 // if tag found, store data in readGroup, return success
\r
673 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
675 // make sure array is more than big enough
\r
676 char newTagData[originalTagDataLength + sizeof(value)];
\r
678 // copy original tag data up til desired tag
\r
679 const unsigned int beginningTagDataLength = numBytesParsed;
\r
680 newTagDataLength += beginningTagDataLength;
\r
681 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
683 // copy new VALUE in place of current tag data
\r
684 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
686 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
\r
688 // skip to next tag (if tag for removal is last, return true)
\r
689 const char* pTagStorageType = pTagData - 1;
\r
690 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
692 // copy everything from current tag (the next one after tag for removal) to end
\r
693 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
694 const unsigned int endTagOffset = beginningTagDataLength + sizeof(float);
\r
695 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
696 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
698 // ensure null-terminator
\r
699 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
701 // save new tag data
\r
702 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
706 // tag not found, attempt AddTag
\r
707 else return AddTag(tag, type, value);
\r
710 // get "NM" tag data - originally contributed by Aaron Quinlan
\r
711 // stores data in 'editDistance', returns success/fail
\r
713 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
\r
714 return GetTag("NM", (uint32_t&)editDistance);
\r
717 // get "RG" tag data
\r
718 // stores data in 'readGroup', returns success/fail
\r
720 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
\r
721 return GetTag("RG", readGroup);
\r
725 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
\r
727 // make sure tag data exists
\r
728 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
731 // localize the tag data
\r
732 char* pTagData = (char*)TagData.data();
\r
733 const unsigned int tagDataLength = TagData.size();
\r
734 unsigned int numBytesParsed = 0;
\r
736 // if tag found, store data in readGroup, return success
\r
737 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
738 const unsigned int dataLength = strlen(pTagData);
\r
739 destination.clear();
\r
740 destination.resize(dataLength);
\r
741 memcpy( (char*)destination.data(), pTagData, dataLength );
\r
745 // tag not found, return failure
\r
750 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
\r
752 // make sure tag data exists
\r
753 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
756 // localize the tag data
\r
757 char* pTagData = (char*)TagData.data();
\r
758 const unsigned int tagDataLength = TagData.size();
\r
759 unsigned int numBytesParsed = 0;
\r
761 // if tag found, determine data byte-length, store data in readGroup, return success
\r
762 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
764 // determine data byte-length
\r
765 const char type = *(pTagData - 1);
\r
766 int destinationLength = 0;
\r
772 destinationLength = 1;
\r
778 destinationLength = 2;
\r
784 destinationLength = 4;
\r
787 // unsupported type for integer destination (float or var-length strings)
\r
791 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
794 // unknown tag type
\r
796 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
800 // store in destination
\r
802 memcpy(&destination, pTagData, destinationLength);
\r
806 // tag not found, return failure
\r
811 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
\r
812 return GetTag(tag, (uint32_t&)destination);
\r
816 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
\r
818 // make sure tag data exists
\r
819 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
822 // localize the tag data
\r
823 char* pTagData = (char*)TagData.data();
\r
824 const unsigned int tagDataLength = TagData.size();
\r
825 unsigned int numBytesParsed = 0;
\r
827 // if tag found, determine data byte-length, store data in readGroup, return success
\r
828 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
829 //pTagData += numBytesParsed;
\r
831 // determine data byte-length
\r
832 const char type = *(pTagData - 1);
\r
833 int destinationLength = 0;
\r
840 destinationLength = 1;
\r
846 destinationLength = 2;
\r
853 destinationLength = 4;
\r
856 // unsupported type (var-length strings)
\r
859 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
862 // unknown tag type
\r
864 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
868 // store in destination
\r
870 memcpy(&destination, pTagData, destinationLength);
\r
874 // tag not found, return failure
\r
879 bool BamAlignment::RemoveTag(const std::string& tag) {
\r
881 // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
\r
882 // also, return false if no data present to remove
\r
883 if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
\r
885 // localize the tag data
\r
886 char* pOriginalTagData = (char*)TagData.data();
\r
887 char* pTagData = pOriginalTagData;
\r
888 const unsigned int originalTagDataLength = TagData.size();
\r
889 unsigned int newTagDataLength = 0;
\r
890 unsigned int numBytesParsed = 0;
\r
892 // if tag found, store data in readGroup, return success
\r
893 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
895 char newTagData[originalTagDataLength];
\r
897 // copy original tag data up til desired tag
\r
899 numBytesParsed -= 3;
\r
900 const unsigned int beginningTagDataLength = numBytesParsed;
\r
901 newTagDataLength += beginningTagDataLength;
\r
902 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
904 // skip to next tag (if tag for removal is last, return true)
\r
905 const char* pTagStorageType = pTagData + 2;
\r
907 numBytesParsed += 3;
\r
908 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
910 // copy everything from current tag (the next one after tag for removal) to end
\r
911 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
912 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
913 memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
\r
915 // save new tag data
\r
916 TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
\r
920 // tag not found, no removal - return failure
\r
925 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {
\r
927 while ( numBytesParsed < tagDataLength ) {
\r
929 const char* pTagType = pTagData;
\r
930 const char* pTagStorageType = pTagData + 2;
\r
932 numBytesParsed += 3;
\r
934 // check the current tag, return true on match
\r
935 if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 )
\r
938 // get the storage class and find the next tag
\r
939 if ( *pTagStorageType == '\0' ) return false;
\r
940 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
\r
941 if ( *pTagData == '\0' ) return false;
\r
944 // checked all tags, none match
\r
949 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
\r
951 switch(storageType) {
\r
962 numBytesParsed += 2;
\r
969 numBytesParsed += 4;
\r
979 // increment for null-terminator
\r
986 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
\r
994 } // namespace BamTools
\r