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_CMATCH = 0;
\r
50 const int BAM_CINS = 1;
\r
51 const int BAM_CDEL = 2;
\r
52 const int BAM_CREF_SKIP = 3;
\r
53 const int BAM_CSOFT_CLIP = 4;
\r
54 const int BAM_CHARD_CLIP = 5;
\r
55 const int BAM_CPAD = 6;
\r
56 const int BAM_CIGAR_SHIFT = 4;
\r
57 const int BAM_CIGAR_MASK = ((1 << BAM_CIGAR_SHIFT) - 1);
\r
58 const int BAM_CORE_SIZE = 32;
\r
59 const int BT_SIZEOF_INT = 4;
\r
63 struct BamAlignment {
\r
65 // constructors & destructor
\r
68 BamAlignment(const BamAlignment& other);
\r
69 ~BamAlignment(void);
\r
71 // Queries against alignment flags
\r
73 bool IsDuplicate(void) const; // Returns true if this read is a PCR duplicate
\r
74 bool IsFailedQC(void) const; // Returns true if this read failed quality control
\r
75 bool IsFirstMate(void) const; // Returns true if alignment is first mate on read
\r
76 bool IsMapped(void) const; // Returns true if alignment is mapped
\r
77 bool IsMateMapped(void) const; // Returns true if alignment's mate is mapped
\r
78 bool IsMateReverseStrand(void) const; // Returns true if alignment's mate mapped to reverse strand
\r
79 bool IsPaired(void) const; // Returns true if alignment part of paired-end read
\r
80 bool IsPrimaryAlignment(void) const; // Returns true if reported position is primary alignment
\r
81 bool IsProperPair(void) const; // Returns true if alignment is part of read that satisfied paired-end resolution
\r
82 bool IsReverseStrand(void) const; // Returns true if alignment mapped to reverse strand
\r
83 bool IsSecondMate(void) const; // Returns true if alignment is second mate on read
\r
85 // Manipulate alignment flags
\r
87 void SetIsDuplicate(bool ok); // Sets "PCR duplicate" flag
\r
88 void SetIsFailedQC(bool ok); // Sets "failed quality control" flag
\r
89 void SetIsFirstMate(bool ok); // Sets "alignment is first mate" flag
\r
90 void SetIsMateUnmapped(bool ok); // Sets "alignment's mate is mapped" flag
\r
91 void SetIsMateReverseStrand(bool ok); // Sets "alignment's mate mapped to reverse strand" flag
\r
92 void SetIsPaired(bool ok); // Sets "alignment part of paired-end read" flag
\r
93 void SetIsProperPair(bool ok); // Sets "alignment is part of read that satisfied paired-end resolution" flag
\r
94 void SetIsReverseStrand(bool ok); // Sets "alignment mapped to reverse strand" flag
\r
95 void SetIsSecondaryAlignment(bool ok); // Sets "position is primary alignment" flag
\r
96 void SetIsSecondMate(bool ok); // Sets "alignment is second mate on read" flag
\r
97 void SetIsUnmapped(bool ok); // Sets "alignment is mapped" flag
\r
99 // Tag data access methods
\r
101 // -------------------------------------------------------------------------------------
\r
102 // N.B. - The following tag-modifying methods may not be used on BamAlignments fetched
\r
103 // using BamReader::GetNextAlignmentCore(). Attempting to use them will not result in
\r
104 // error message (to keep output clean) but will ALWAYS return false. Only user-
\r
105 // generated BamAlignments or those retrieved using BamReader::GetNextAlignment() are valid.
\r
107 // add tag data (create new TAG entry with TYPE and VALUE)
\r
108 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
109 // returns true if new data added, false if error or TAG already exists
\r
110 // N.B. - will NOT modify existing tag. Use EditTag() instead
\r
111 bool AddTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
112 bool AddTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
113 bool AddTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
114 bool AddTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
116 // edit tag data (sets existing TAG with TYPE to VALUE or adds new TAG if not already present)
\r
117 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
118 // returns true if edit was successfaul, false if error
\r
119 bool EditTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
120 bool EditTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
121 bool EditTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
122 bool EditTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
124 // specific tag data access methods - these only remain for legacy support
\r
125 bool GetEditDistance(uint32_t& editDistance) const; // get "NM" tag data (implemented as GetTag("NM", editDistance))
\r
126 bool GetReadGroup(std::string& readGroup) const; // get "RG" tag data (implemented as GetTag("RG", readGroup))
\r
128 // generic tag data access methods
\r
129 bool GetTag(const std::string& tag, std::string& destination) const; // access variable-length char or hex strings
\r
130 bool GetTag(const std::string& tag, uint32_t& destination) const; // access unsigned integer data
\r
131 bool GetTag(const std::string& tag, int32_t& destination) const; // access signed integer data
\r
132 bool GetTag(const std::string& tag, float& destination) const; // access floating point data
\r
135 // returns true if removal was successful, false if error
\r
136 // N.B. - returns false if TAG does not exist (no removal can occur)
\r
137 bool RemoveTag(const std::string& tag);
\r
139 // Additional data access methods
\r
141 // calculates alignment end position, based on starting position and CIGAR operations
\r
142 // @zeroBased - if true, returns 0-based coordinate; else returns 1-based
\r
143 int GetEndPosition(bool usePadded = false, bool zeroBased = true) const;
\r
145 // 'internal' utility methods
\r
147 static bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed);
\r
148 static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);
\r
152 std::string Name; // Read name
\r
153 int32_t Length; // Query length
\r
154 std::string QueryBases; // 'Original' sequence (as reported from sequencing machine)
\r
155 std::string AlignedBases; // 'Aligned' sequence (includes any indels, padding, clipping)
\r
156 std::string Qualities; // FASTQ qualities (ASCII characters, not numeric values)
\r
157 std::string TagData; // Tag data (accessor methods will pull the requested information out)
\r
158 int32_t RefID; // ID number for reference sequence
\r
159 int32_t Position; // Position (0-based) where alignment starts
\r
160 uint16_t Bin; // Bin in BAM file where this alignment resides
\r
161 uint16_t MapQuality; // Mapping quality score
\r
162 uint32_t AlignmentFlag; // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate
\r
163 std::vector<CigarOp> CigarData; // CIGAR operations for this alignment
\r
164 int32_t MateRefID; // ID number for reference sequence where alignment's mate was aligned
\r
165 int32_t MatePosition; // Position (0-based) where alignment's mate starts
\r
166 int32_t InsertSize; // Mate-pair insert size
\r
170 struct BamAlignmentSupportData {
\r
173 std::string AllCharData;
\r
174 uint32_t BlockLength;
\r
175 uint32_t NumCigarOperations;
\r
176 uint32_t QueryNameLength;
\r
177 uint32_t QuerySequenceLength;
\r
181 BamAlignmentSupportData(void)
\r
183 , NumCigarOperations(0)
\r
184 , QueryNameLength(0)
\r
185 , QuerySequenceLength(0)
\r
186 , HasCoreOnly(false)
\r
190 // contains raw character data & lengths
\r
191 BamAlignmentSupportData SupportData;
\r
193 // allow these classes access to BamAlignment private members (SupportData)
\r
194 // but client code should not need to touch this data
\r
195 friend class BamReader;
\r
196 friend class BamWriter;
\r
198 // Alignment flag query constants
\r
199 // Use the get/set methods above instead
\r
204 , MATE_UNMAPPED = 8
\r
206 , MATE_REVERSE = 32
\r
211 , DUPLICATE = 1024
\r
215 // ----------------------------------------------------------------
\r
216 // Auxiliary data structs & typedefs
\r
221 char Type; // Operation type (MIDNSHP)
\r
222 uint32_t Length; // Operation length (number of bases)
\r
225 CigarOp(const char type = '\0',
\r
226 const uint32_t length = 0)
\r
235 std::string RefName; // Name of reference sequence
\r
236 int32_t RefLength; // Length of reference sequence
\r
237 bool RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence
\r
240 RefData(const int32_t& length = 0,
\r
242 : RefLength(length)
\r
243 , RefHasAlignments(ok)
\r
247 typedef std::vector<RefData> RefVector;
\r
248 typedef std::vector<BamAlignment> BamAlignmentVector;
\r
259 BamRegion(const int& leftID = -1,
\r
260 const int& leftPos = -1,
\r
261 const int& rightID = -1,
\r
262 const int& rightPos = -1)
\r
263 : LeftRefID(leftID)
\r
264 , LeftPosition(leftPos)
\r
265 , RightRefID(rightID)
\r
266 , RightPosition(rightPos)
\r
269 // member functions
\r
270 void clear(void) { LeftRefID = -1; LeftPosition = -1; RightRefID = -1; RightPosition = -1; }
\r
271 bool isLeftBoundSpecified(void) const { return ( LeftRefID != -1 && LeftPosition != -1 ); }
\r
272 bool isNull(void) const { return ( !isLeftBoundSpecified() && !isRightBoundSpecified() ); }
\r
273 bool isRightBoundSpecified(void) const { return ( RightRefID != -1 && RightPosition != -1 ); }
\r
276 // ----------------------------------------------------------------
\r
277 // Added: 3-35-2010 DWB
\r
278 // Fixed: Routines to provide endian-correctness
\r
279 // ----------------------------------------------------------------
\r
281 // returns true if system is big endian
\r
282 inline bool SystemIsBigEndian(void) {
\r
283 const uint16_t one = 0x0001;
\r
284 return ((*(char*) &one) == 0 );
\r
287 // swaps endianness of 16-bit value 'in place'
\r
288 inline void SwapEndian_16(int16_t& x) {
\r
289 x = ((x >> 8) | (x << 8));
\r
292 inline void SwapEndian_16(uint16_t& x) {
\r
293 x = ((x >> 8) | (x << 8));
\r
296 // swaps endianness of 32-bit value 'in-place'
\r
297 inline void SwapEndian_32(int32_t& x) {
\r
299 ((x << 8) & 0x00FF0000) |
\r
300 ((x >> 8) & 0x0000FF00) |
\r
305 inline void SwapEndian_32(uint32_t& x) {
\r
307 ((x << 8) & 0x00FF0000) |
\r
308 ((x >> 8) & 0x0000FF00) |
\r
313 // swaps endianness of 64-bit value 'in-place'
\r
314 inline void SwapEndian_64(int64_t& x) {
\r
316 ((x << 40) & 0x00FF000000000000ll) |
\r
317 ((x << 24) & 0x0000FF0000000000ll) |
\r
318 ((x << 8) & 0x000000FF00000000ll) |
\r
319 ((x >> 8) & 0x00000000FF000000ll) |
\r
320 ((x >> 24) & 0x0000000000FF0000ll) |
\r
321 ((x >> 40) & 0x000000000000FF00ll) |
\r
326 inline void SwapEndian_64(uint64_t& x) {
\r
328 ((x << 40) & 0x00FF000000000000ll) |
\r
329 ((x << 24) & 0x0000FF0000000000ll) |
\r
330 ((x << 8) & 0x000000FF00000000ll) |
\r
331 ((x >> 8) & 0x00000000FF000000ll) |
\r
332 ((x >> 24) & 0x0000000000FF0000ll) |
\r
333 ((x >> 40) & 0x000000000000FF00ll) |
\r
338 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)
\r
339 inline void SwapEndian_16p(char* data) {
\r
340 uint16_t& value = (uint16_t&)*data;
\r
341 SwapEndian_16(value);
\r
344 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)
\r
345 inline void SwapEndian_32p(char* data) {
\r
346 uint32_t& value = (uint32_t&)*data;
\r
347 SwapEndian_32(value);
\r
350 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)
\r
351 inline void SwapEndian_64p(char* data) {
\r
352 uint64_t& value = (uint64_t&)*data;
\r
353 SwapEndian_64(value);
\r
356 inline bool FileExists(const std::string& filename) {
\r
357 std::ifstream f(filename.c_str(), std::ifstream::in);
\r
361 // ----------------------------------------------------------------
\r
362 // BamAlignment member methods
\r
364 // constructors & destructor
\r
365 inline BamAlignment::BamAlignment(void) { }
\r
367 inline BamAlignment::BamAlignment(const BamAlignment& other)
\r
369 , Length(other.Length)
\r
370 , QueryBases(other.QueryBases)
\r
371 , AlignedBases(other.AlignedBases)
\r
372 , Qualities(other.Qualities)
\r
373 , TagData(other.TagData)
\r
374 , RefID(other.RefID)
\r
375 , Position(other.Position)
\r
377 , MapQuality(other.MapQuality)
\r
378 , AlignmentFlag(other.AlignmentFlag)
\r
379 , CigarData(other.CigarData)
\r
380 , MateRefID(other.MateRefID)
\r
381 , MatePosition(other.MatePosition)
\r
382 , InsertSize(other.InsertSize)
\r
383 , SupportData(other.SupportData)
\r
386 inline BamAlignment::~BamAlignment(void) { }
\r
388 // Queries against alignment flags
\r
389 inline bool BamAlignment::IsDuplicate(void) const { return ( (AlignmentFlag & DUPLICATE) != 0 ); }
\r
390 inline bool BamAlignment::IsFailedQC(void) const { return ( (AlignmentFlag & QC_FAILED) != 0 ); }
\r
391 inline bool BamAlignment::IsFirstMate(void) const { return ( (AlignmentFlag & READ_1) != 0 ); }
\r
392 inline bool BamAlignment::IsMapped(void) const { return ( (AlignmentFlag & UNMAPPED) == 0 ); }
\r
393 inline bool BamAlignment::IsMateMapped(void) const { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }
\r
394 inline bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE) != 0 ); }
\r
395 inline bool BamAlignment::IsPaired(void) const { return ( (AlignmentFlag & PAIRED) != 0 ); }
\r
396 inline bool BamAlignment::IsPrimaryAlignment(void) const { return ( (AlignmentFlag & SECONDARY) == 0 ); }
\r
397 inline bool BamAlignment::IsProperPair(void) const { return ( (AlignmentFlag & PROPER_PAIR) != 0 ); }
\r
398 inline bool BamAlignment::IsReverseStrand(void) const { return ( (AlignmentFlag & REVERSE) != 0 ); }
\r
399 inline bool BamAlignment::IsSecondMate(void) const { return ( (AlignmentFlag & READ_2) != 0 ); }
\r
401 // Manipulate alignment flags
\r
402 inline void BamAlignment::SetIsDuplicate(bool ok) { if (ok) AlignmentFlag |= DUPLICATE; else AlignmentFlag &= ~DUPLICATE; }
\r
403 inline void BamAlignment::SetIsFailedQC(bool ok) { if (ok) AlignmentFlag |= QC_FAILED; else AlignmentFlag &= ~QC_FAILED; }
\r
404 inline void BamAlignment::SetIsFirstMate(bool ok) { if (ok) AlignmentFlag |= READ_1; else AlignmentFlag &= ~READ_1; }
\r
405 inline void BamAlignment::SetIsMateUnmapped(bool ok) { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }
\r
406 inline void BamAlignment::SetIsMateReverseStrand(bool ok) { if (ok) AlignmentFlag |= MATE_REVERSE; else AlignmentFlag &= ~MATE_REVERSE; }
\r
407 inline void BamAlignment::SetIsPaired(bool ok) { if (ok) AlignmentFlag |= PAIRED; else AlignmentFlag &= ~PAIRED; }
\r
408 inline void BamAlignment::SetIsProperPair(bool ok) { if (ok) AlignmentFlag |= PROPER_PAIR; else AlignmentFlag &= ~PROPER_PAIR; }
\r
409 inline void BamAlignment::SetIsReverseStrand(bool ok) { if (ok) AlignmentFlag |= REVERSE; else AlignmentFlag &= ~REVERSE; }
\r
410 inline void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY; else AlignmentFlag &= ~SECONDARY; }
\r
411 inline void BamAlignment::SetIsSecondMate(bool ok) { if (ok) AlignmentFlag |= READ_2; else AlignmentFlag &= ~READ_2; }
\r
412 inline void BamAlignment::SetIsUnmapped(bool ok) { if (ok) AlignmentFlag |= UNMAPPED; else AlignmentFlag &= ~UNMAPPED; }
\r
414 // calculates alignment end position, based on starting position and CIGAR operations
\r
416 int BamAlignment::GetEndPosition(bool usePadded, bool zeroBased) const {
\r
418 // initialize alignment end to starting position
\r
419 int alignEnd = Position;
\r
421 // iterate over cigar operations
\r
422 std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
\r
423 std::vector<CigarOp>::const_iterator cigarEnd = CigarData.end();
\r
424 for ( ; cigarIter != cigarEnd; ++cigarIter) {
\r
425 const char cigarType = (*cigarIter).Type;
\r
426 if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {
\r
427 alignEnd += (*cigarIter).Length;
\r
429 else if ( usePadded && cigarType == 'I' ) {
\r
430 alignEnd += (*cigarIter).Length;
\r
434 // adjust for zeroBased, if necessary
\r
436 return alignEnd - 1;
\r
442 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
444 if ( SupportData.HasCoreOnly ) return false;
\r
445 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
446 if ( type != "Z" && type != "H" ) return false;
\r
448 // localize the tag data
\r
449 char* pTagData = (char*)TagData.data();
\r
450 const unsigned int tagDataLength = TagData.size();
\r
451 unsigned int numBytesParsed = 0;
\r
453 // if tag already exists, return false
\r
454 // use EditTag explicitly instead
\r
455 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
457 // otherwise, copy tag data to temp buffer
\r
458 std::string newTag = tag + type + value;
\r
459 const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
\r
460 char originalTagData[newTagDataLength];
\r
461 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
464 strcat(originalTagData + tagDataLength, newTag.data()); // removes original null-term, appends newTag + null-term
\r
466 // store temp buffer back in TagData
\r
467 const char* newTagData = (const char*)originalTagData;
\r
468 TagData.assign(newTagData, newTagDataLength);
\r
475 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
477 if ( SupportData.HasCoreOnly ) return false;
\r
478 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
479 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
481 // localize the tag data
\r
482 char* pTagData = (char*)TagData.data();
\r
483 const unsigned int tagDataLength = TagData.size();
\r
484 unsigned int numBytesParsed = 0;
\r
486 // if tag already exists, return false
\r
487 // use EditTag explicitly instead
\r
488 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
490 // otherwise, convert value to string
\r
491 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
494 // copy original tag data to temp buffer
\r
495 std::string newTag = tag + type;
\r
496 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
\r
497 char originalTagData[newTagDataLength];
\r
498 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
501 strcat(originalTagData + tagDataLength, newTag.data());
\r
502 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
\r
504 // store temp buffer back in TagData
\r
505 const char* newTagData = (const char*)originalTagData;
\r
506 TagData.assign(newTagData, newTagDataLength);
\r
513 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
514 return AddTag(tag, type, (const uint32_t&)value);
\r
518 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
\r
520 if ( SupportData.HasCoreOnly ) return false;
\r
521 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
522 if ( type == "Z" || type == "H" ) return false;
\r
524 // localize the tag data
\r
525 char* pTagData = (char*)TagData.data();
\r
526 const unsigned int tagDataLength = TagData.size();
\r
527 unsigned int numBytesParsed = 0;
\r
529 // if tag already exists, return false
\r
530 // use EditTag explicitly instead
\r
531 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
533 // otherwise, convert value to string
\r
534 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
537 // copy original tag data to temp buffer
\r
538 std::string newTag = tag + type;
\r
539 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
\r
540 char originalTagData[newTagDataLength];
\r
541 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
544 strcat(originalTagData + tagDataLength, newTag.data());
\r
545 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
\r
547 // store temp buffer back in TagData
\r
548 const char* newTagData = (const char*)originalTagData;
\r
549 TagData.assign(newTagData, newTagDataLength);
\r
556 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
558 if ( SupportData.HasCoreOnly ) return false;
\r
559 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
560 if ( type != "Z" && type != "H" ) return false;
\r
562 // localize the tag data
\r
563 char* pOriginalTagData = (char*)TagData.data();
\r
564 char* pTagData = pOriginalTagData;
\r
565 const unsigned int originalTagDataLength = TagData.size();
\r
567 unsigned int newTagDataLength = 0;
\r
568 unsigned int numBytesParsed = 0;
\r
570 // if tag found, store data in readGroup, return success
\r
571 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
573 // make sure array is more than big enough
\r
574 char newTagData[originalTagDataLength + value.size()];
\r
576 // copy original tag data up til desired tag
\r
577 const unsigned int beginningTagDataLength = numBytesParsed;
\r
578 newTagDataLength += beginningTagDataLength;
\r
579 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
581 // copy new VALUE in place of current tag data
\r
582 const unsigned int dataLength = strlen(value.c_str());
\r
583 memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
\r
585 // skip to next tag (if tag for removal is last, return true)
\r
586 const char* pTagStorageType = pTagData - 1;
\r
587 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
589 // copy everything from current tag (the next one after tag for removal) to end
\r
590 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
591 const unsigned int endTagOffset = beginningTagDataLength + dataLength + 1;
\r
592 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
593 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
595 // ensure null-terminator
\r
596 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
598 // save new tag data
\r
599 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
603 // tag not found, attempt AddTag
\r
604 else return AddTag(tag, type, value);
\r
608 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
610 if ( SupportData.HasCoreOnly ) return false;
\r
611 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
612 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
614 // localize the tag data
\r
615 char* pOriginalTagData = (char*)TagData.data();
\r
616 char* pTagData = pOriginalTagData;
\r
617 const unsigned int originalTagDataLength = TagData.size();
\r
619 unsigned int newTagDataLength = 0;
\r
620 unsigned int numBytesParsed = 0;
\r
622 // if tag found, store data in readGroup, return success
\r
623 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
625 // make sure array is more than big enough
\r
626 char newTagData[originalTagDataLength + sizeof(value)];
\r
628 // copy original tag data up til desired tag
\r
629 const unsigned int beginningTagDataLength = numBytesParsed;
\r
630 newTagDataLength += beginningTagDataLength;
\r
631 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
633 // copy new VALUE in place of current tag data
\r
634 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
636 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
\r
638 // skip to next tag (if tag for removal is last, return true)
\r
639 const char* pTagStorageType = pTagData - 1;
\r
640 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
642 // copy everything from current tag (the next one after tag for removal) to end
\r
643 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
644 const unsigned int endTagOffset = beginningTagDataLength + sizeof(unsigned int);
\r
645 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
646 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
648 // ensure null-terminator
\r
649 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
651 // save new tag data
\r
652 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
656 // tag not found, attempt AddTag
\r
657 else return AddTag(tag, type, value);
\r
661 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
662 return EditTag(tag, type, (const uint32_t&)value);
\r
666 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
\r
668 if ( SupportData.HasCoreOnly ) return false;
\r
669 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
670 if ( type == "Z" || type == "H" ) return false;
\r
672 // localize the tag data
\r
673 char* pOriginalTagData = (char*)TagData.data();
\r
674 char* pTagData = pOriginalTagData;
\r
675 const unsigned int originalTagDataLength = TagData.size();
\r
677 unsigned int newTagDataLength = 0;
\r
678 unsigned int numBytesParsed = 0;
\r
680 // if tag found, store data in readGroup, return success
\r
681 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
683 // make sure array is more than big enough
\r
684 char newTagData[originalTagDataLength + sizeof(value)];
\r
686 // copy original tag data up til desired tag
\r
687 const unsigned int beginningTagDataLength = numBytesParsed;
\r
688 newTagDataLength += beginningTagDataLength;
\r
689 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
691 // copy new VALUE in place of current tag data
\r
692 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
694 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
\r
696 // skip to next tag (if tag for removal is last, return true)
\r
697 const char* pTagStorageType = pTagData - 1;
\r
698 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
700 // copy everything from current tag (the next one after tag for removal) to end
\r
701 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
702 const unsigned int endTagOffset = beginningTagDataLength + sizeof(float);
\r
703 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
704 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
706 // ensure null-terminator
\r
707 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
709 // save new tag data
\r
710 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
714 // tag not found, attempt AddTag
\r
715 else return AddTag(tag, type, value);
\r
718 // get "NM" tag data - originally contributed by Aaron Quinlan
\r
719 // stores data in 'editDistance', returns success/fail
\r
721 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
\r
722 return GetTag("NM", (uint32_t&)editDistance);
\r
725 // get "RG" tag data
\r
726 // stores data in 'readGroup', returns success/fail
\r
728 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
\r
729 return GetTag("RG", readGroup);
\r
733 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
\r
735 // make sure tag data exists
\r
736 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
739 // localize the tag data
\r
740 char* pTagData = (char*)TagData.data();
\r
741 const unsigned int tagDataLength = TagData.size();
\r
742 unsigned int numBytesParsed = 0;
\r
744 // if tag found, store data in readGroup, return success
\r
745 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
746 const unsigned int dataLength = strlen(pTagData);
\r
747 destination.clear();
\r
748 destination.resize(dataLength);
\r
749 memcpy( (char*)destination.data(), pTagData, dataLength );
\r
753 // tag not found, return failure
\r
758 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
\r
760 // make sure tag data exists
\r
761 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
764 // localize the tag data
\r
765 char* pTagData = (char*)TagData.data();
\r
766 const unsigned int tagDataLength = TagData.size();
\r
767 unsigned int numBytesParsed = 0;
\r
769 // if tag found, determine data byte-length, store data in readGroup, return success
\r
770 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
772 // determine data byte-length
\r
773 const char type = *(pTagData - 1);
\r
774 int destinationLength = 0;
\r
780 destinationLength = 1;
\r
786 destinationLength = 2;
\r
792 destinationLength = 4;
\r
795 // unsupported type for integer destination (float or var-length strings)
\r
799 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
802 // unknown tag type
\r
804 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
808 // store in destination
\r
810 memcpy(&destination, pTagData, destinationLength);
\r
814 // tag not found, return failure
\r
819 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
\r
820 return GetTag(tag, (uint32_t&)destination);
\r
824 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
\r
826 // make sure tag data exists
\r
827 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
830 // localize the tag data
\r
831 char* pTagData = (char*)TagData.data();
\r
832 const unsigned int tagDataLength = TagData.size();
\r
833 unsigned int numBytesParsed = 0;
\r
835 // if tag found, determine data byte-length, store data in readGroup, return success
\r
836 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
837 //pTagData += numBytesParsed;
\r
839 // determine data byte-length
\r
840 const char type = *(pTagData - 1);
\r
841 int destinationLength = 0;
\r
848 destinationLength = 1;
\r
854 destinationLength = 2;
\r
861 destinationLength = 4;
\r
864 // unsupported type (var-length strings)
\r
867 fprintf(stderr, "ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
870 // unknown tag type
\r
872 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
876 // store in destination
\r
878 memcpy(&destination, pTagData, destinationLength);
\r
882 // tag not found, return failure
\r
887 bool BamAlignment::RemoveTag(const std::string& tag) {
\r
889 // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
\r
890 // also, return false if no data present to remove
\r
891 if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
\r
893 // localize the tag data
\r
894 char* pOriginalTagData = (char*)TagData.data();
\r
895 char* pTagData = pOriginalTagData;
\r
896 const unsigned int originalTagDataLength = TagData.size();
\r
897 unsigned int newTagDataLength = 0;
\r
898 unsigned int numBytesParsed = 0;
\r
900 // if tag found, store data in readGroup, return success
\r
901 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
903 char newTagData[originalTagDataLength];
\r
905 // copy original tag data up til desired tag
\r
907 numBytesParsed -= 3;
\r
908 const unsigned int beginningTagDataLength = numBytesParsed;
\r
909 newTagDataLength += beginningTagDataLength;
\r
910 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
912 // skip to next tag (if tag for removal is last, return true)
\r
913 const char* pTagStorageType = pTagData + 2;
\r
915 numBytesParsed += 3;
\r
916 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
918 // copy everything from current tag (the next one after tag for removal) to end
\r
919 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
920 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
921 memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
\r
923 // save new tag data
\r
924 TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
\r
928 // tag not found, no removal - return failure
\r
933 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {
\r
935 while ( numBytesParsed < tagDataLength ) {
\r
937 const char* pTagType = pTagData;
\r
938 const char* pTagStorageType = pTagData + 2;
\r
940 numBytesParsed += 3;
\r
942 // check the current tag, return true on match
\r
943 if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 )
\r
946 // get the storage class and find the next tag
\r
947 if ( *pTagStorageType == '\0' ) return false;
\r
948 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
\r
949 if ( *pTagData == '\0' ) return false;
\r
952 // checked all tags, none match
\r
957 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
\r
959 switch(storageType) {
\r
970 numBytesParsed += 2;
\r
977 numBytesParsed += 4;
\r
987 // increment for null-terminator
\r
994 fprintf(stderr, "ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
\r
1002 } // namespace BamTools
\r
1004 #endif // BAMAUX_H
\r