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: 27 July 2010 (DB)
\r
7 // ---------------------------------------------------------------------------
\r
8 // Provides the basic constants, data structures, etc. for using BAM files
\r
9 // ***************************************************************************
\r
21 #include <exception>
\r
27 // Platform-specific type definitions
\r
28 #ifndef BAMTOOLS_TYPES
\r
29 #define BAMTOOLS_TYPES
\r
31 typedef char int8_t;
\r
32 typedef unsigned char uint8_t;
\r
33 typedef short int16_t;
\r
34 typedef unsigned short uint16_t;
\r
35 typedef int int32_t;
\r
36 typedef unsigned int uint32_t;
\r
37 typedef long long int64_t;
\r
38 typedef unsigned long long uint64_t;
\r
42 #endif // BAMTOOLS_TYPES
\r
44 namespace BamTools {
\r
47 const int BAM_CORE_SIZE = 32;
\r
48 const int BAM_CMATCH = 0;
\r
49 const int BAM_CINS = 1;
\r
50 const int BAM_CDEL = 2;
\r
51 const int BAM_CREF_SKIP = 3;
\r
52 const int BAM_CSOFT_CLIP = 4;
\r
53 const int BAM_CHARD_CLIP = 5;
\r
54 const int BAM_CPAD = 6;
\r
55 const int BAM_CIGAR_SHIFT = 4;
\r
56 const int BAM_CIGAR_MASK = ((1 << BAM_CIGAR_SHIFT) - 1);
\r
58 // BAM index constants
\r
59 const int MAX_BIN = 37450; // =(8^6-1)/7+1
\r
60 const int BAM_MIN_CHUNK_GAP = 32768;
\r
61 const int BAM_LIDX_SHIFT = 14;
\r
63 // Explicit variable sizes
\r
64 const int BT_SIZEOF_INT = 4;
\r
68 struct BamAlignment {
\r
70 // constructors & destructor
\r
73 BamAlignment(const BamAlignment& other);
\r
74 ~BamAlignment(void);
\r
76 // Queries against alignment flags
\r
78 bool IsDuplicate(void) const; // Returns true if this read is a PCR duplicate
\r
79 bool IsFailedQC(void) const; // Returns true if this read failed quality control
\r
80 bool IsFirstMate(void) const; // Returns true if alignment is first mate on read
\r
81 bool IsMapped(void) const; // Returns true if alignment is mapped
\r
82 bool IsMateMapped(void) const; // Returns true if alignment's mate is mapped
\r
83 bool IsMateReverseStrand(void) const; // Returns true if alignment's mate mapped to reverse strand
\r
84 bool IsPaired(void) const; // Returns true if alignment part of paired-end read
\r
85 bool IsPrimaryAlignment(void) const; // Returns true if reported position is primary alignment
\r
86 bool IsProperPair(void) const; // Returns true if alignment is part of read that satisfied paired-end resolution
\r
87 bool IsReverseStrand(void) const; // Returns true if alignment mapped to reverse strand
\r
88 bool IsSecondMate(void) const; // Returns true if alignment is second mate on read
\r
90 // Manipulate alignment flags
\r
92 void SetIsDuplicate(bool ok); // Sets "PCR duplicate" flag
\r
93 void SetIsFailedQC(bool ok); // Sets "failed quality control" flag
\r
94 void SetIsFirstMate(bool ok); // Sets "alignment is first mate" flag
\r
95 void SetIsMateUnmapped(bool ok); // Sets "alignment's mate is mapped" flag
\r
96 void SetIsMateReverseStrand(bool ok); // Sets "alignment's mate mapped to reverse strand" flag
\r
97 void SetIsPaired(bool ok); // Sets "alignment part of paired-end read" flag
\r
98 void SetIsProperPair(bool ok); // Sets "alignment is part of read that satisfied paired-end resolution" flag
\r
99 void SetIsReverseStrand(bool ok); // Sets "alignment mapped to reverse strand" flag
\r
100 void SetIsSecondaryAlignment(bool ok); // Sets "position is primary alignment" flag
\r
101 void SetIsSecondMate(bool ok); // Sets "alignment is second mate on read" flag
\r
102 void SetIsUnmapped(bool ok); // Sets "alignment is mapped" flag
\r
104 // Tag data access methods
\r
106 // -------------------------------------------------------------------------------------
\r
107 // N.B. - The following tag-modifying methods may not be used on BamAlignments fetched
\r
108 // using BamReader::GetNextAlignmentCore(). Attempting to use them will not result in
\r
109 // error message (to keep output clean) but will ALWAYS return false. Only user-
\r
110 // generated BamAlignments or those retrieved using BamReader::GetNextAlignment() are valid.
\r
112 // add tag data (create new TAG entry with TYPE and VALUE)
\r
113 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
114 // returns true if new data added, false if error or TAG already exists
\r
115 // N.B. - will NOT modify existing tag. Use EditTag() instead
\r
116 bool AddTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
117 bool AddTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
118 bool AddTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
119 bool AddTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
121 // edit tag data (sets existing TAG with TYPE to VALUE or adds new TAG if not already present)
\r
122 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
123 // returns true if edit was successfaul, false if error
\r
124 bool EditTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
125 bool EditTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
126 bool EditTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
127 bool EditTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
129 // specific tag data access methods - these only remain for legacy support
\r
130 bool GetEditDistance(uint32_t& editDistance) const; // get "NM" tag data (implemented as GetTag("NM", editDistance))
\r
131 bool GetReadGroup(std::string& readGroup) const; // get "RG" tag data (implemented as GetTag("RG", readGroup))
\r
133 // generic tag data access methods
\r
134 bool GetTag(const std::string& tag, std::string& destination) const; // access variable-length char or hex strings
\r
135 bool GetTag(const std::string& tag, uint32_t& destination) const; // access unsigned integer data
\r
136 bool GetTag(const std::string& tag, int32_t& destination) const; // access signed integer data
\r
137 bool GetTag(const std::string& tag, float& destination) const; // access floating point data
\r
140 // returns true if removal was successful, false if error
\r
141 // N.B. - returns false if TAG does not exist (no removal can occur)
\r
142 bool RemoveTag(const std::string& tag);
\r
144 // Additional data access methods
\r
146 int GetEndPosition(bool usePadded = false) const; // calculates alignment end position, based on starting position and CIGAR operations
\r
148 // 'internal' utility methods
\r
150 static bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed);
\r
151 static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);
\r
155 std::string Name; // Read name
\r
156 int32_t Length; // Query length
\r
157 std::string QueryBases; // 'Original' sequence (as reported from sequencing machine)
\r
158 std::string AlignedBases; // 'Aligned' sequence (includes any indels, padding, clipping)
\r
159 std::string Qualities; // FASTQ qualities (ASCII characters, not numeric values)
\r
160 std::string TagData; // Tag data (accessor methods will pull the requested information out)
\r
161 int32_t RefID; // ID number for reference sequence
\r
162 int32_t Position; // Position (0-based) where alignment starts
\r
163 uint16_t Bin; // Bin in BAM file where this alignment resides
\r
164 uint16_t MapQuality; // Mapping quality score
\r
165 uint32_t AlignmentFlag; // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate
\r
166 std::vector<CigarOp> CigarData; // CIGAR operations for this alignment
\r
167 int32_t MateRefID; // ID number for reference sequence where alignment's mate was aligned
\r
168 int32_t MatePosition; // Position (0-based) where alignment's mate starts
\r
169 int32_t InsertSize; // Mate-pair insert size
\r
173 struct BamAlignmentSupportData {
\r
176 std::string AllCharData;
\r
177 uint32_t BlockLength;
\r
178 uint32_t NumCigarOperations;
\r
179 uint32_t QueryNameLength;
\r
180 uint32_t QuerySequenceLength;
\r
184 BamAlignmentSupportData(void)
\r
186 , NumCigarOperations(0)
\r
187 , QueryNameLength(0)
\r
188 , QuerySequenceLength(0)
\r
189 , HasCoreOnly(false)
\r
193 // contains raw character data & lengths
\r
194 BamAlignmentSupportData SupportData;
\r
196 // allow these classes access to BamAlignment private members (SupportData)
\r
197 // but client code should not need to touch this data
\r
198 friend class BamReader;
\r
199 friend class BamWriter;
\r
201 // Alignment flag query constants
\r
202 // Use the get/set methods above instead
\r
207 , MATE_UNMAPPED = 8
\r
209 , MATE_REVERSE = 32
\r
214 , DUPLICATE = 1024
\r
218 // ----------------------------------------------------------------
\r
219 // Auxiliary data structs & typedefs
\r
224 char Type; // Operation type (MIDNSHP)
\r
225 uint32_t Length; // Operation length (number of bases)
\r
228 CigarOp(const char type = '\0',
\r
229 const uint32_t length = 0)
\r
238 std::string RefName; // Name of reference sequence
\r
239 int32_t RefLength; // Length of reference sequence
\r
240 bool RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence
\r
243 RefData(const int32_t& length = 0,
\r
245 : RefLength(length)
\r
246 , RefHasAlignments(ok)
\r
250 typedef std::vector<RefData> RefVector;
\r
251 typedef std::vector<BamAlignment> BamAlignmentVector;
\r
262 BamRegion(const int& leftID = -1,
\r
263 const int& leftPos = -1,
\r
264 const int& rightID = -1,
\r
265 const int& rightPos = -1)
\r
266 : LeftRefID(leftID)
\r
267 , LeftPosition(leftPos)
\r
268 , RightRefID(rightID)
\r
269 , RightPosition(rightPos)
\r
273 // ----------------------------------------------------------------
\r
274 // Added: 3-35-2010 DWB
\r
275 // Fixed: Routines to provide endian-correctness
\r
276 // ----------------------------------------------------------------
\r
278 // returns true if system is big endian
\r
279 inline bool SystemIsBigEndian(void) {
\r
280 const uint16_t one = 0x0001;
\r
281 return ((*(char*) &one) == 0 );
\r
284 // swaps endianness of 16-bit value 'in place'
\r
285 inline void SwapEndian_16(int16_t& x) {
\r
286 x = ((x >> 8) | (x << 8));
\r
289 inline void SwapEndian_16(uint16_t& x) {
\r
290 x = ((x >> 8) | (x << 8));
\r
293 // swaps endianness of 32-bit value 'in-place'
\r
294 inline void SwapEndian_32(int32_t& x) {
\r
296 ((x << 8) & 0x00FF0000) |
\r
297 ((x >> 8) & 0x0000FF00) |
\r
302 inline void SwapEndian_32(uint32_t& x) {
\r
304 ((x << 8) & 0x00FF0000) |
\r
305 ((x >> 8) & 0x0000FF00) |
\r
310 // swaps endianness of 64-bit value 'in-place'
\r
311 inline void SwapEndian_64(int64_t& x) {
\r
313 ((x << 40) & 0x00FF000000000000ll) |
\r
314 ((x << 24) & 0x0000FF0000000000ll) |
\r
315 ((x << 8) & 0x000000FF00000000ll) |
\r
316 ((x >> 8) & 0x00000000FF000000ll) |
\r
317 ((x >> 24) & 0x0000000000FF0000ll) |
\r
318 ((x >> 40) & 0x000000000000FF00ll) |
\r
323 inline void SwapEndian_64(uint64_t& x) {
\r
325 ((x << 40) & 0x00FF000000000000ll) |
\r
326 ((x << 24) & 0x0000FF0000000000ll) |
\r
327 ((x << 8) & 0x000000FF00000000ll) |
\r
328 ((x >> 8) & 0x00000000FF000000ll) |
\r
329 ((x >> 24) & 0x0000000000FF0000ll) |
\r
330 ((x >> 40) & 0x000000000000FF00ll) |
\r
335 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)
\r
336 inline void SwapEndian_16p(char* data) {
\r
337 uint16_t& value = (uint16_t&)*data;
\r
338 SwapEndian_16(value);
\r
341 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)
\r
342 inline void SwapEndian_32p(char* data) {
\r
343 uint32_t& value = (uint32_t&)*data;
\r
344 SwapEndian_32(value);
\r
347 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)
\r
348 inline void SwapEndian_64p(char* data) {
\r
349 uint64_t& value = (uint64_t&)*data;
\r
350 SwapEndian_64(value);
\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) 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
429 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
431 if ( SupportData.HasCoreOnly ) return false;
\r
432 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
433 if ( type != "Z" && type != "H" ) return false;
\r
435 // localize the tag data
\r
436 char* pTagData = (char*)TagData.data();
\r
437 const unsigned int tagDataLength = TagData.size();
\r
438 unsigned int numBytesParsed = 0;
\r
440 // if tag already exists, return false
\r
441 // use EditTag explicitly instead
\r
442 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
444 // otherwise, copy tag data to temp buffer
\r
445 std::string newTag = tag + type + value;
\r
446 const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
\r
447 char originalTagData[newTagDataLength];
\r
448 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
451 strcat(originalTagData + tagDataLength, newTag.data()); // removes original null-term, appends newTag + null-term
\r
453 // store temp buffer back in TagData
\r
454 const char* newTagData = (const char*)originalTagData;
\r
455 TagData.assign(newTagData, newTagDataLength);
\r
462 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
464 if ( SupportData.HasCoreOnly ) return false;
\r
465 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
466 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
468 // localize the tag data
\r
469 char* pTagData = (char*)TagData.data();
\r
470 const unsigned int tagDataLength = TagData.size();
\r
471 unsigned int numBytesParsed = 0;
\r
473 // if tag already exists, return false
\r
474 // use EditTag explicitly instead
\r
475 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
477 // otherwise, convert value to string
\r
478 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
481 // copy original tag data to temp buffer
\r
482 std::string newTag = tag + type;
\r
483 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
\r
484 char originalTagData[newTagDataLength];
\r
485 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
488 strcat(originalTagData + tagDataLength, newTag.data());
\r
489 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
\r
491 // store temp buffer back in TagData
\r
492 const char* newTagData = (const char*)originalTagData;
\r
493 TagData.assign(newTagData, newTagDataLength);
\r
500 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
501 return AddTag(tag, type, (const uint32_t&)value);
\r
505 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
\r
507 if ( SupportData.HasCoreOnly ) return false;
\r
508 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
509 if ( type == "Z" || type == "H" ) return false;
\r
511 // localize the tag data
\r
512 char* pTagData = (char*)TagData.data();
\r
513 const unsigned int tagDataLength = TagData.size();
\r
514 unsigned int numBytesParsed = 0;
\r
516 // if tag already exists, return false
\r
517 // use EditTag explicitly instead
\r
518 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
520 // otherwise, convert value to string
\r
521 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
524 // copy original tag data to temp buffer
\r
525 std::string newTag = tag + type;
\r
526 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
\r
527 char originalTagData[newTagDataLength];
\r
528 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
531 strcat(originalTagData + tagDataLength, newTag.data());
\r
532 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
\r
534 // store temp buffer back in TagData
\r
535 const char* newTagData = (const char*)originalTagData;
\r
536 TagData.assign(newTagData, newTagDataLength);
\r
543 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
545 if ( SupportData.HasCoreOnly ) return false;
\r
546 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
547 if ( type != "Z" && type != "H" ) return false;
\r
549 // localize the tag data
\r
550 char* pOriginalTagData = (char*)TagData.data();
\r
551 char* pTagData = pOriginalTagData;
\r
552 const unsigned int originalTagDataLength = TagData.size();
\r
554 unsigned int newTagDataLength = 0;
\r
555 unsigned int numBytesParsed = 0;
\r
557 // if tag found, store data in readGroup, return success
\r
558 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
560 // make sure array is more than big enough
\r
561 char newTagData[originalTagDataLength + value.size()];
\r
563 // copy original tag data up til desired tag
\r
564 const unsigned int beginningTagDataLength = numBytesParsed;
\r
565 newTagDataLength += beginningTagDataLength;
\r
566 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
568 // copy new VALUE in place of current tag data
\r
569 const unsigned int dataLength = strlen(value.c_str());
\r
570 memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
\r
572 // skip to next tag (if tag for removal is last, return true)
\r
573 const char* pTagStorageType = pTagData - 1;
\r
574 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
576 // copy everything from current tag (the next one after tag for removal) to end
\r
577 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
578 const unsigned int endTagOffset = beginningTagDataLength + dataLength + 1;
\r
579 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
580 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
582 // ensure null-terminator
\r
583 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
585 // save new tag data
\r
586 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
590 // tag not found, attempt AddTag
\r
591 else return AddTag(tag, type, value);
\r
595 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
597 if ( SupportData.HasCoreOnly ) return false;
\r
598 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
599 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
601 // localize the tag data
\r
602 char* pOriginalTagData = (char*)TagData.data();
\r
603 char* pTagData = pOriginalTagData;
\r
604 const unsigned int originalTagDataLength = TagData.size();
\r
606 unsigned int newTagDataLength = 0;
\r
607 unsigned int numBytesParsed = 0;
\r
609 // if tag found, store data in readGroup, return success
\r
610 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
612 // make sure array is more than big enough
\r
613 char newTagData[originalTagDataLength + sizeof(value)];
\r
615 // copy original tag data up til desired tag
\r
616 const unsigned int beginningTagDataLength = numBytesParsed;
\r
617 newTagDataLength += beginningTagDataLength;
\r
618 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
620 // copy new VALUE in place of current tag data
\r
621 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
623 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
\r
625 // skip to next tag (if tag for removal is last, return true)
\r
626 const char* pTagStorageType = pTagData - 1;
\r
627 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
629 // copy everything from current tag (the next one after tag for removal) to end
\r
630 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
631 const unsigned int endTagOffset = beginningTagDataLength + sizeof(unsigned int);
\r
632 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
633 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
635 // ensure null-terminator
\r
636 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
638 // save new tag data
\r
639 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
643 // tag not found, attempt AddTag
\r
644 else return AddTag(tag, type, value);
\r
648 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
649 return EditTag(tag, type, (const uint32_t&)value);
\r
653 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
\r
655 if ( SupportData.HasCoreOnly ) return false;
\r
656 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
657 if ( type == "Z" || type == "H" ) return false;
\r
659 // localize the tag data
\r
660 char* pOriginalTagData = (char*)TagData.data();
\r
661 char* pTagData = pOriginalTagData;
\r
662 const unsigned int originalTagDataLength = TagData.size();
\r
664 unsigned int newTagDataLength = 0;
\r
665 unsigned int numBytesParsed = 0;
\r
667 // if tag found, store data in readGroup, return success
\r
668 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
670 // make sure array is more than big enough
\r
671 char newTagData[originalTagDataLength + sizeof(value)];
\r
673 // copy original tag data up til desired tag
\r
674 const unsigned int beginningTagDataLength = numBytesParsed;
\r
675 newTagDataLength += beginningTagDataLength;
\r
676 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
678 // copy new VALUE in place of current tag data
\r
679 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
681 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
\r
683 // skip to next tag (if tag for removal is last, return true)
\r
684 const char* pTagStorageType = pTagData - 1;
\r
685 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
687 // copy everything from current tag (the next one after tag for removal) to end
\r
688 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
689 const unsigned int endTagOffset = beginningTagDataLength + sizeof(float);
\r
690 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
691 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
693 // ensure null-terminator
\r
694 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
696 // save new tag data
\r
697 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
701 // tag not found, attempt AddTag
\r
702 else return AddTag(tag, type, value);
\r
705 // get "NM" tag data - originally contributed by Aaron Quinlan
\r
706 // stores data in 'editDistance', returns success/fail
\r
708 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
\r
709 return GetTag("NM", (uint32_t&)editDistance);
\r
712 // get "RG" tag data
\r
713 // stores data in 'readGroup', returns success/fail
\r
715 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
\r
716 return GetTag("RG", readGroup);
\r
720 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
\r
722 // make sure tag data exists
\r
723 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
726 // localize the tag data
\r
727 char* pTagData = (char*)TagData.data();
\r
728 const unsigned int tagDataLength = TagData.size();
\r
729 unsigned int numBytesParsed = 0;
\r
731 // if tag found, store data in readGroup, return success
\r
732 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
733 const unsigned int dataLength = strlen(pTagData);
\r
734 destination.clear();
\r
735 destination.resize(dataLength);
\r
736 memcpy( (char*)destination.data(), pTagData, dataLength );
\r
740 // tag not found, return failure
\r
745 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
\r
747 // make sure tag data exists
\r
748 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
751 // localize the tag data
\r
752 char* pTagData = (char*)TagData.data();
\r
753 const unsigned int tagDataLength = TagData.size();
\r
754 unsigned int numBytesParsed = 0;
\r
756 // if tag found, determine data byte-length, store data in readGroup, return success
\r
757 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
759 // determine data byte-length
\r
760 const char type = *(pTagData - 1);
\r
761 int destinationLength = 0;
\r
767 destinationLength = 1;
\r
773 destinationLength = 2;
\r
779 destinationLength = 4;
\r
782 // unsupported type for integer destination (float or var-length strings)
\r
786 printf("ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
789 // unknown tag type
\r
791 printf("ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
795 // store in destination
\r
797 memcpy(&destination, pTagData, destinationLength);
\r
801 // tag not found, return failure
\r
806 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
\r
807 return GetTag(tag, (uint32_t&)destination);
\r
811 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
\r
813 // make sure tag data exists
\r
814 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
817 // localize the tag data
\r
818 char* pTagData = (char*)TagData.data();
\r
819 const unsigned int tagDataLength = TagData.size();
\r
820 unsigned int numBytesParsed = 0;
\r
822 // if tag found, determine data byte-length, store data in readGroup, return success
\r
823 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
824 //pTagData += numBytesParsed;
\r
826 // determine data byte-length
\r
827 const char type = *(pTagData - 1);
\r
828 int destinationLength = 0;
\r
835 destinationLength = 1;
\r
841 destinationLength = 2;
\r
848 destinationLength = 4;
\r
851 // unsupported type (var-length strings)
\r
854 printf("ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
857 // unknown tag type
\r
859 printf("ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
863 // store in destination
\r
865 memcpy(&destination, pTagData, destinationLength);
\r
869 // tag not found, return failure
\r
874 bool BamAlignment::RemoveTag(const std::string& tag) {
\r
876 // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
\r
877 // also, return false if no data present to remove
\r
878 if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
\r
880 // localize the tag data
\r
881 char* pOriginalTagData = (char*)TagData.data();
\r
882 char* pTagData = pOriginalTagData;
\r
883 const unsigned int originalTagDataLength = TagData.size();
\r
884 unsigned int newTagDataLength = 0;
\r
885 unsigned int numBytesParsed = 0;
\r
887 // if tag found, store data in readGroup, return success
\r
888 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
890 char newTagData[originalTagDataLength];
\r
892 // copy original tag data up til desired tag
\r
894 numBytesParsed -= 3;
\r
895 const unsigned int beginningTagDataLength = numBytesParsed;
\r
896 newTagDataLength += beginningTagDataLength;
\r
897 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
899 // skip to next tag (if tag for removal is last, return true)
\r
900 const char* pTagStorageType = pTagData + 2;
\r
902 numBytesParsed += 3;
\r
903 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
905 // copy everything from current tag (the next one after tag for removal) to end
\r
906 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
907 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
908 memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
\r
910 // save new tag data
\r
911 TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
\r
915 // tag not found, no removal - return failure
\r
920 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {
\r
922 while ( numBytesParsed < tagDataLength ) {
\r
924 const char* pTagType = pTagData;
\r
925 const char* pTagStorageType = pTagData + 2;
\r
927 numBytesParsed += 3;
\r
929 // check the current tag, return true on match
\r
930 if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 )
\r
933 // get the storage class and find the next tag
\r
934 if ( *pTagStorageType == '\0' ) return false;
\r
935 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
\r
936 if ( *pTagData == '\0' ) return false;
\r
939 // checked all tags, none match
\r
944 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
\r
946 switch(storageType) {
\r
957 numBytesParsed += 2;
\r
964 numBytesParsed += 4;
\r
974 // increment for null-terminator
\r
981 printf("ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
\r
989 } // namespace BamTools
\r