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: 3 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_CORE_SIZE = 32;
\r
50 const int BAM_CMATCH = 0;
\r
51 const int BAM_CINS = 1;
\r
52 const int BAM_CDEL = 2;
\r
53 const int BAM_CREF_SKIP = 3;
\r
54 const int BAM_CSOFT_CLIP = 4;
\r
55 const int BAM_CHARD_CLIP = 5;
\r
56 const int BAM_CPAD = 6;
\r
57 const int BAM_CIGAR_SHIFT = 4;
\r
58 const int BAM_CIGAR_MASK = ((1 << BAM_CIGAR_SHIFT) - 1);
\r
60 // BAM index constants
\r
61 const int MAX_BIN = 37450; // =(8^6-1)/7+1
\r
62 const int BAM_MIN_CHUNK_GAP = 32768;
\r
63 const int BAM_LIDX_SHIFT = 14;
\r
65 // Explicit variable sizes
\r
66 const int BT_SIZEOF_INT = 4;
\r
70 struct BamAlignment {
\r
72 // constructors & destructor
\r
75 BamAlignment(const BamAlignment& other);
\r
76 ~BamAlignment(void);
\r
78 // Queries against alignment flags
\r
80 bool IsDuplicate(void) const; // Returns true if this read is a PCR duplicate
\r
81 bool IsFailedQC(void) const; // Returns true if this read failed quality control
\r
82 bool IsFirstMate(void) const; // Returns true if alignment is first mate on read
\r
83 bool IsMapped(void) const; // Returns true if alignment is mapped
\r
84 bool IsMateMapped(void) const; // Returns true if alignment's mate is mapped
\r
85 bool IsMateReverseStrand(void) const; // Returns true if alignment's mate mapped to reverse strand
\r
86 bool IsPaired(void) const; // Returns true if alignment part of paired-end read
\r
87 bool IsPrimaryAlignment(void) const; // Returns true if reported position is primary alignment
\r
88 bool IsProperPair(void) const; // Returns true if alignment is part of read that satisfied paired-end resolution
\r
89 bool IsReverseStrand(void) const; // Returns true if alignment mapped to reverse strand
\r
90 bool IsSecondMate(void) const; // Returns true if alignment is second mate on read
\r
92 // Manipulate alignment flags
\r
94 void SetIsDuplicate(bool ok); // Sets "PCR duplicate" flag
\r
95 void SetIsFailedQC(bool ok); // Sets "failed quality control" flag
\r
96 void SetIsFirstMate(bool ok); // Sets "alignment is first mate" flag
\r
97 void SetIsMateUnmapped(bool ok); // Sets "alignment's mate is mapped" flag
\r
98 void SetIsMateReverseStrand(bool ok); // Sets "alignment's mate mapped to reverse strand" flag
\r
99 void SetIsPaired(bool ok); // Sets "alignment part of paired-end read" flag
\r
100 void SetIsProperPair(bool ok); // Sets "alignment is part of read that satisfied paired-end resolution" flag
\r
101 void SetIsReverseStrand(bool ok); // Sets "alignment mapped to reverse strand" flag
\r
102 void SetIsSecondaryAlignment(bool ok); // Sets "position is primary alignment" flag
\r
103 void SetIsSecondMate(bool ok); // Sets "alignment is second mate on read" flag
\r
104 void SetIsUnmapped(bool ok); // Sets "alignment is mapped" flag
\r
106 // Tag data access methods
\r
108 // -------------------------------------------------------------------------------------
\r
109 // N.B. - The following tag-modifying methods may not be used on BamAlignments fetched
\r
110 // using BamReader::GetNextAlignmentCore(). Attempting to use them will not result in
\r
111 // error message (to keep output clean) but will ALWAYS return false. Only user-
\r
112 // generated BamAlignments or those retrieved using BamReader::GetNextAlignment() are valid.
\r
114 // add tag data (create new TAG entry with TYPE and VALUE)
\r
115 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
116 // returns true if new data added, false if error or TAG already exists
\r
117 // N.B. - will NOT modify existing tag. Use EditTag() instead
\r
118 bool AddTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
119 bool AddTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
120 bool AddTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
121 bool AddTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
123 // edit tag data (sets existing TAG with TYPE to VALUE or adds new TAG if not already present)
\r
124 // TYPE is one of {A, i, f, Z, H} depending on VALUE - see SAM/BAM spec for details
\r
125 // returns true if edit was successfaul, false if error
\r
126 bool EditTag(const std::string& tag, const std::string& type, const std::string& value); // type must be Z or H
\r
127 bool EditTag(const std::string& tag, const std::string& type, const uint32_t& value); // type must be A or i
\r
128 bool EditTag(const std::string& tag, const std::string& type, const int32_t& value); // type must be A or i
\r
129 bool EditTag(const std::string& tag, const std::string& type, const float& value); // type must be A, i, or f
\r
131 // specific tag data access methods - these only remain for legacy support
\r
132 bool GetEditDistance(uint32_t& editDistance) const; // get "NM" tag data (implemented as GetTag("NM", editDistance))
\r
133 bool GetReadGroup(std::string& readGroup) const; // get "RG" tag data (implemented as GetTag("RG", readGroup))
\r
135 // generic tag data access methods
\r
136 bool GetTag(const std::string& tag, std::string& destination) const; // access variable-length char or hex strings
\r
137 bool GetTag(const std::string& tag, uint32_t& destination) const; // access unsigned integer data
\r
138 bool GetTag(const std::string& tag, int32_t& destination) const; // access signed integer data
\r
139 bool GetTag(const std::string& tag, float& destination) const; // access floating point data
\r
142 // returns true if removal was successful, false if error
\r
143 // N.B. - returns false if TAG does not exist (no removal can occur)
\r
144 bool RemoveTag(const std::string& tag);
\r
146 // Additional data access methods
\r
148 int GetEndPosition(bool usePadded = false) const; // calculates alignment end position, based on starting position and CIGAR operations
\r
150 // 'internal' utility methods
\r
152 static bool FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed);
\r
153 static bool SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed);
\r
157 std::string Name; // Read name
\r
158 int32_t Length; // Query length
\r
159 std::string QueryBases; // 'Original' sequence (as reported from sequencing machine)
\r
160 std::string AlignedBases; // 'Aligned' sequence (includes any indels, padding, clipping)
\r
161 std::string Qualities; // FASTQ qualities (ASCII characters, not numeric values)
\r
162 std::string TagData; // Tag data (accessor methods will pull the requested information out)
\r
163 int32_t RefID; // ID number for reference sequence
\r
164 int32_t Position; // Position (0-based) where alignment starts
\r
165 uint16_t Bin; // Bin in BAM file where this alignment resides
\r
166 uint16_t MapQuality; // Mapping quality score
\r
167 uint32_t AlignmentFlag; // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate
\r
168 std::vector<CigarOp> CigarData; // CIGAR operations for this alignment
\r
169 int32_t MateRefID; // ID number for reference sequence where alignment's mate was aligned
\r
170 int32_t MatePosition; // Position (0-based) where alignment's mate starts
\r
171 int32_t InsertSize; // Mate-pair insert size
\r
175 struct BamAlignmentSupportData {
\r
178 std::string AllCharData;
\r
179 uint32_t BlockLength;
\r
180 uint32_t NumCigarOperations;
\r
181 uint32_t QueryNameLength;
\r
182 uint32_t QuerySequenceLength;
\r
186 BamAlignmentSupportData(void)
\r
188 , NumCigarOperations(0)
\r
189 , QueryNameLength(0)
\r
190 , QuerySequenceLength(0)
\r
191 , HasCoreOnly(false)
\r
195 // contains raw character data & lengths
\r
196 BamAlignmentSupportData SupportData;
\r
198 // allow these classes access to BamAlignment private members (SupportData)
\r
199 // but client code should not need to touch this data
\r
200 friend class BamReader;
\r
201 friend class BamWriter;
\r
203 // Alignment flag query constants
\r
204 // Use the get/set methods above instead
\r
209 , MATE_UNMAPPED = 8
\r
211 , MATE_REVERSE = 32
\r
216 , DUPLICATE = 1024
\r
220 // ----------------------------------------------------------------
\r
221 // Auxiliary data structs & typedefs
\r
226 char Type; // Operation type (MIDNSHP)
\r
227 uint32_t Length; // Operation length (number of bases)
\r
230 CigarOp(const char type = '\0',
\r
231 const uint32_t length = 0)
\r
240 std::string RefName; // Name of reference sequence
\r
241 int32_t RefLength; // Length of reference sequence
\r
242 bool RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence
\r
245 RefData(const int32_t& length = 0,
\r
247 : RefLength(length)
\r
248 , RefHasAlignments(ok)
\r
252 typedef std::vector<RefData> RefVector;
\r
253 typedef std::vector<BamAlignment> BamAlignmentVector;
\r
264 BamRegion(const int& leftID = -1,
\r
265 const int& leftPos = -1,
\r
266 const int& rightID = -1,
\r
267 const int& rightPos = -1)
\r
268 : LeftRefID(leftID)
\r
269 , LeftPosition(leftPos)
\r
270 , RightRefID(rightID)
\r
271 , RightPosition(rightPos)
\r
275 // ----------------------------------------------------------------
\r
276 // Added: 3-35-2010 DWB
\r
277 // Fixed: Routines to provide endian-correctness
\r
278 // ----------------------------------------------------------------
\r
280 // returns true if system is big endian
\r
281 inline bool SystemIsBigEndian(void) {
\r
282 const uint16_t one = 0x0001;
\r
283 return ((*(char*) &one) == 0 );
\r
286 // swaps endianness of 16-bit value 'in place'
\r
287 inline void SwapEndian_16(int16_t& x) {
\r
288 x = ((x >> 8) | (x << 8));
\r
291 inline void SwapEndian_16(uint16_t& x) {
\r
292 x = ((x >> 8) | (x << 8));
\r
295 // swaps endianness of 32-bit value 'in-place'
\r
296 inline void SwapEndian_32(int32_t& x) {
\r
298 ((x << 8) & 0x00FF0000) |
\r
299 ((x >> 8) & 0x0000FF00) |
\r
304 inline void SwapEndian_32(uint32_t& x) {
\r
306 ((x << 8) & 0x00FF0000) |
\r
307 ((x >> 8) & 0x0000FF00) |
\r
312 // swaps endianness of 64-bit value 'in-place'
\r
313 inline void SwapEndian_64(int64_t& x) {
\r
315 ((x << 40) & 0x00FF000000000000ll) |
\r
316 ((x << 24) & 0x0000FF0000000000ll) |
\r
317 ((x << 8) & 0x000000FF00000000ll) |
\r
318 ((x >> 8) & 0x00000000FF000000ll) |
\r
319 ((x >> 24) & 0x0000000000FF0000ll) |
\r
320 ((x >> 40) & 0x000000000000FF00ll) |
\r
325 inline void SwapEndian_64(uint64_t& x) {
\r
327 ((x << 40) & 0x00FF000000000000ll) |
\r
328 ((x << 24) & 0x0000FF0000000000ll) |
\r
329 ((x << 8) & 0x000000FF00000000ll) |
\r
330 ((x >> 8) & 0x00000000FF000000ll) |
\r
331 ((x >> 24) & 0x0000000000FF0000ll) |
\r
332 ((x >> 40) & 0x000000000000FF00ll) |
\r
337 // swaps endianness of 'next 2 bytes' in a char buffer (in-place)
\r
338 inline void SwapEndian_16p(char* data) {
\r
339 uint16_t& value = (uint16_t&)*data;
\r
340 SwapEndian_16(value);
\r
343 // swaps endianness of 'next 4 bytes' in a char buffer (in-place)
\r
344 inline void SwapEndian_32p(char* data) {
\r
345 uint32_t& value = (uint32_t&)*data;
\r
346 SwapEndian_32(value);
\r
349 // swaps endianness of 'next 8 bytes' in a char buffer (in-place)
\r
350 inline void SwapEndian_64p(char* data) {
\r
351 uint64_t& value = (uint64_t&)*data;
\r
352 SwapEndian_64(value);
\r
355 inline bool FileExists(const std::string& filename) {
\r
356 std::ifstream f(filename.c_str(), std::ifstream::in);
\r
360 // ----------------------------------------------------------------
\r
361 // BamAlignment member methods
\r
363 // constructors & destructor
\r
364 inline BamAlignment::BamAlignment(void) { }
\r
366 inline BamAlignment::BamAlignment(const BamAlignment& other)
\r
368 , Length(other.Length)
\r
369 , QueryBases(other.QueryBases)
\r
370 , AlignedBases(other.AlignedBases)
\r
371 , Qualities(other.Qualities)
\r
372 , TagData(other.TagData)
\r
373 , RefID(other.RefID)
\r
374 , Position(other.Position)
\r
376 , MapQuality(other.MapQuality)
\r
377 , AlignmentFlag(other.AlignmentFlag)
\r
378 , CigarData(other.CigarData)
\r
379 , MateRefID(other.MateRefID)
\r
380 , MatePosition(other.MatePosition)
\r
381 , InsertSize(other.InsertSize)
\r
382 , SupportData(other.SupportData)
\r
385 inline BamAlignment::~BamAlignment(void) { }
\r
387 // Queries against alignment flags
\r
388 inline bool BamAlignment::IsDuplicate(void) const { return ( (AlignmentFlag & DUPLICATE) != 0 ); }
\r
389 inline bool BamAlignment::IsFailedQC(void) const { return ( (AlignmentFlag & QC_FAILED) != 0 ); }
\r
390 inline bool BamAlignment::IsFirstMate(void) const { return ( (AlignmentFlag & READ_1) != 0 ); }
\r
391 inline bool BamAlignment::IsMapped(void) const { return ( (AlignmentFlag & UNMAPPED) == 0 ); }
\r
392 inline bool BamAlignment::IsMateMapped(void) const { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }
\r
393 inline bool BamAlignment::IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE) != 0 ); }
\r
394 inline bool BamAlignment::IsPaired(void) const { return ( (AlignmentFlag & PAIRED) != 0 ); }
\r
395 inline bool BamAlignment::IsPrimaryAlignment(void) const { return ( (AlignmentFlag & SECONDARY) == 0 ); }
\r
396 inline bool BamAlignment::IsProperPair(void) const { return ( (AlignmentFlag & PROPER_PAIR) != 0 ); }
\r
397 inline bool BamAlignment::IsReverseStrand(void) const { return ( (AlignmentFlag & REVERSE) != 0 ); }
\r
398 inline bool BamAlignment::IsSecondMate(void) const { return ( (AlignmentFlag & READ_2) != 0 ); }
\r
400 // Manipulate alignment flags
\r
401 inline void BamAlignment::SetIsDuplicate(bool ok) { if (ok) AlignmentFlag |= DUPLICATE; else AlignmentFlag &= ~DUPLICATE; }
\r
402 inline void BamAlignment::SetIsFailedQC(bool ok) { if (ok) AlignmentFlag |= QC_FAILED; else AlignmentFlag &= ~QC_FAILED; }
\r
403 inline void BamAlignment::SetIsFirstMate(bool ok) { if (ok) AlignmentFlag |= READ_1; else AlignmentFlag &= ~READ_1; }
\r
404 inline void BamAlignment::SetIsMateUnmapped(bool ok) { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }
\r
405 inline void BamAlignment::SetIsMateReverseStrand(bool ok) { if (ok) AlignmentFlag |= MATE_REVERSE; else AlignmentFlag &= ~MATE_REVERSE; }
\r
406 inline void BamAlignment::SetIsPaired(bool ok) { if (ok) AlignmentFlag |= PAIRED; else AlignmentFlag &= ~PAIRED; }
\r
407 inline void BamAlignment::SetIsProperPair(bool ok) { if (ok) AlignmentFlag |= PROPER_PAIR; else AlignmentFlag &= ~PROPER_PAIR; }
\r
408 inline void BamAlignment::SetIsReverseStrand(bool ok) { if (ok) AlignmentFlag |= REVERSE; else AlignmentFlag &= ~REVERSE; }
\r
409 inline void BamAlignment::SetIsSecondaryAlignment(bool ok) { if (ok) AlignmentFlag |= SECONDARY; else AlignmentFlag &= ~SECONDARY; }
\r
410 inline void BamAlignment::SetIsSecondMate(bool ok) { if (ok) AlignmentFlag |= READ_2; else AlignmentFlag &= ~READ_2; }
\r
411 inline void BamAlignment::SetIsUnmapped(bool ok) { if (ok) AlignmentFlag |= UNMAPPED; else AlignmentFlag &= ~UNMAPPED; }
\r
413 // calculates alignment end position, based on starting position and CIGAR operations
\r
415 int BamAlignment::GetEndPosition(bool usePadded) const {
\r
417 // initialize alignment end to starting position
\r
418 int alignEnd = Position;
\r
420 // iterate over cigar operations
\r
421 std::vector<CigarOp>::const_iterator cigarIter = CigarData.begin();
\r
422 std::vector<CigarOp>::const_iterator cigarEnd = CigarData.end();
\r
423 for ( ; cigarIter != cigarEnd; ++cigarIter) {
\r
424 const char cigarType = (*cigarIter).Type;
\r
425 if ( cigarType == 'M' || cigarType == 'D' || cigarType == 'N' ) {
\r
426 alignEnd += (*cigarIter).Length;
\r
428 else if ( usePadded && cigarType == 'I' ) {
\r
429 alignEnd += (*cigarIter).Length;
\r
436 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
438 if ( SupportData.HasCoreOnly ) return false;
\r
439 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
440 if ( type != "Z" && type != "H" ) return false;
\r
442 // localize the tag data
\r
443 char* pTagData = (char*)TagData.data();
\r
444 const unsigned int tagDataLength = TagData.size();
\r
445 unsigned int numBytesParsed = 0;
\r
447 // if tag already exists, return false
\r
448 // use EditTag explicitly instead
\r
449 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
451 // otherwise, copy tag data to temp buffer
\r
452 std::string newTag = tag + type + value;
\r
453 const int newTagDataLength = tagDataLength + newTag.size() + 1; // leave room for null-term
\r
454 char originalTagData[newTagDataLength];
\r
455 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
458 strcat(originalTagData + tagDataLength, newTag.data()); // removes original null-term, appends newTag + null-term
\r
460 // store temp buffer back in TagData
\r
461 const char* newTagData = (const char*)originalTagData;
\r
462 TagData.assign(newTagData, newTagDataLength);
\r
469 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
471 if ( SupportData.HasCoreOnly ) return false;
\r
472 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
473 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
475 // localize the tag data
\r
476 char* pTagData = (char*)TagData.data();
\r
477 const unsigned int tagDataLength = TagData.size();
\r
478 unsigned int numBytesParsed = 0;
\r
480 // if tag already exists, return false
\r
481 // use EditTag explicitly instead
\r
482 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
484 // otherwise, convert value to string
\r
485 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
488 // copy original tag data to temp buffer
\r
489 std::string newTag = tag + type;
\r
490 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new integer
\r
491 char originalTagData[newTagDataLength];
\r
492 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
495 strcat(originalTagData + tagDataLength, newTag.data());
\r
496 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(unsigned int));
\r
498 // store temp buffer back in TagData
\r
499 const char* newTagData = (const char*)originalTagData;
\r
500 TagData.assign(newTagData, newTagDataLength);
\r
507 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
508 return AddTag(tag, type, (const uint32_t&)value);
\r
512 bool BamAlignment::AddTag(const std::string& tag, const std::string& type, const float& value) {
\r
514 if ( SupportData.HasCoreOnly ) return false;
\r
515 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
516 if ( type == "Z" || type == "H" ) return false;
\r
518 // localize the tag data
\r
519 char* pTagData = (char*)TagData.data();
\r
520 const unsigned int tagDataLength = TagData.size();
\r
521 unsigned int numBytesParsed = 0;
\r
523 // if tag already exists, return false
\r
524 // use EditTag explicitly instead
\r
525 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) return false;
\r
527 // otherwise, convert value to string
\r
528 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
531 // copy original tag data to temp buffer
\r
532 std::string newTag = tag + type;
\r
533 const int newTagDataLength = tagDataLength + newTag.size() + 4; // leave room for new float
\r
534 char originalTagData[newTagDataLength];
\r
535 memcpy(originalTagData, TagData.c_str(), tagDataLength + 1); // '+1' for TagData null-term
\r
538 strcat(originalTagData + tagDataLength, newTag.data());
\r
539 memcpy(originalTagData + tagDataLength + newTag.size(), un.valueBuffer, sizeof(float));
\r
541 // store temp buffer back in TagData
\r
542 const char* newTagData = (const char*)originalTagData;
\r
543 TagData.assign(newTagData, newTagDataLength);
\r
550 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const std::string& value) {
\r
552 if ( SupportData.HasCoreOnly ) return false;
\r
553 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
554 if ( type != "Z" && type != "H" ) return false;
\r
556 // localize the tag data
\r
557 char* pOriginalTagData = (char*)TagData.data();
\r
558 char* pTagData = pOriginalTagData;
\r
559 const unsigned int originalTagDataLength = TagData.size();
\r
561 unsigned int newTagDataLength = 0;
\r
562 unsigned int numBytesParsed = 0;
\r
564 // if tag found, store data in readGroup, return success
\r
565 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
567 // make sure array is more than big enough
\r
568 char newTagData[originalTagDataLength + value.size()];
\r
570 // copy original tag data up til desired tag
\r
571 const unsigned int beginningTagDataLength = numBytesParsed;
\r
572 newTagDataLength += beginningTagDataLength;
\r
573 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
575 // copy new VALUE in place of current tag data
\r
576 const unsigned int dataLength = strlen(value.c_str());
\r
577 memcpy(newTagData + beginningTagDataLength, (char*)value.c_str(), dataLength+1 );
\r
579 // skip to next tag (if tag for removal is last, return true)
\r
580 const char* pTagStorageType = pTagData - 1;
\r
581 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
583 // copy everything from current tag (the next one after tag for removal) to end
\r
584 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
585 const unsigned int endTagOffset = beginningTagDataLength + dataLength + 1;
\r
586 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
587 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
589 // ensure null-terminator
\r
590 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
592 // save new tag data
\r
593 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
597 // tag not found, attempt AddTag
\r
598 else return AddTag(tag, type, value);
\r
602 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const uint32_t& value) {
\r
604 if ( SupportData.HasCoreOnly ) return false;
\r
605 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
606 if ( type == "f" || type == "Z" || type == "H" ) return false;
\r
608 // localize the tag data
\r
609 char* pOriginalTagData = (char*)TagData.data();
\r
610 char* pTagData = pOriginalTagData;
\r
611 const unsigned int originalTagDataLength = TagData.size();
\r
613 unsigned int newTagDataLength = 0;
\r
614 unsigned int numBytesParsed = 0;
\r
616 // if tag found, store data in readGroup, return success
\r
617 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
619 // make sure array is more than big enough
\r
620 char newTagData[originalTagDataLength + sizeof(value)];
\r
622 // copy original tag data up til desired tag
\r
623 const unsigned int beginningTagDataLength = numBytesParsed;
\r
624 newTagDataLength += beginningTagDataLength;
\r
625 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
627 // copy new VALUE in place of current tag data
\r
628 union { unsigned int value; char valueBuffer[sizeof(unsigned int)]; } un;
\r
630 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(unsigned int));
\r
632 // skip to next tag (if tag for removal is last, return true)
\r
633 const char* pTagStorageType = pTagData - 1;
\r
634 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
636 // copy everything from current tag (the next one after tag for removal) to end
\r
637 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
638 const unsigned int endTagOffset = beginningTagDataLength + sizeof(unsigned int);
\r
639 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
640 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
642 // ensure null-terminator
\r
643 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
645 // save new tag data
\r
646 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
650 // tag not found, attempt AddTag
\r
651 else return AddTag(tag, type, value);
\r
655 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const int32_t& value) {
\r
656 return EditTag(tag, type, (const uint32_t&)value);
\r
660 bool BamAlignment::EditTag(const std::string& tag, const std::string& type, const float& value) {
\r
662 if ( SupportData.HasCoreOnly ) return false;
\r
663 if ( tag.size() != 2 || type.size() != 1 ) return false;
\r
664 if ( type == "Z" || type == "H" ) return false;
\r
666 // localize the tag data
\r
667 char* pOriginalTagData = (char*)TagData.data();
\r
668 char* pTagData = pOriginalTagData;
\r
669 const unsigned int originalTagDataLength = TagData.size();
\r
671 unsigned int newTagDataLength = 0;
\r
672 unsigned int numBytesParsed = 0;
\r
674 // if tag found, store data in readGroup, return success
\r
675 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
677 // make sure array is more than big enough
\r
678 char newTagData[originalTagDataLength + sizeof(value)];
\r
680 // copy original tag data up til desired tag
\r
681 const unsigned int beginningTagDataLength = numBytesParsed;
\r
682 newTagDataLength += beginningTagDataLength;
\r
683 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
685 // copy new VALUE in place of current tag data
\r
686 union { float value; char valueBuffer[sizeof(float)]; } un;
\r
688 memcpy(newTagData + beginningTagDataLength, un.valueBuffer, sizeof(float));
\r
690 // skip to next tag (if tag for removal is last, return true)
\r
691 const char* pTagStorageType = pTagData - 1;
\r
692 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
694 // copy everything from current tag (the next one after tag for removal) to end
\r
695 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
696 const unsigned int endTagOffset = beginningTagDataLength + sizeof(float);
\r
697 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
698 memcpy(newTagData + endTagOffset, pTagData, endTagDataLength);
\r
700 // ensure null-terminator
\r
701 newTagData[ endTagOffset + endTagDataLength + 1 ] = 0;
\r
703 // save new tag data
\r
704 TagData.assign(newTagData, endTagOffset + endTagDataLength);
\r
708 // tag not found, attempt AddTag
\r
709 else return AddTag(tag, type, value);
\r
712 // get "NM" tag data - originally contributed by Aaron Quinlan
\r
713 // stores data in 'editDistance', returns success/fail
\r
715 bool BamAlignment::GetEditDistance(uint32_t& editDistance) const {
\r
716 return GetTag("NM", (uint32_t&)editDistance);
\r
719 // get "RG" tag data
\r
720 // stores data in 'readGroup', returns success/fail
\r
722 bool BamAlignment::GetReadGroup(std::string& readGroup) const {
\r
723 return GetTag("RG", readGroup);
\r
727 bool BamAlignment::GetTag(const std::string& tag, std::string& destination) const {
\r
729 // make sure tag data exists
\r
730 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
733 // localize the tag data
\r
734 char* pTagData = (char*)TagData.data();
\r
735 const unsigned int tagDataLength = TagData.size();
\r
736 unsigned int numBytesParsed = 0;
\r
738 // if tag found, store data in readGroup, return success
\r
739 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
740 const unsigned int dataLength = strlen(pTagData);
\r
741 destination.clear();
\r
742 destination.resize(dataLength);
\r
743 memcpy( (char*)destination.data(), pTagData, dataLength );
\r
747 // tag not found, return failure
\r
752 bool BamAlignment::GetTag(const std::string& tag, uint32_t& destination) const {
\r
754 // make sure tag data exists
\r
755 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
758 // localize the tag data
\r
759 char* pTagData = (char*)TagData.data();
\r
760 const unsigned int tagDataLength = TagData.size();
\r
761 unsigned int numBytesParsed = 0;
\r
763 // if tag found, determine data byte-length, store data in readGroup, return success
\r
764 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
766 // determine data byte-length
\r
767 const char type = *(pTagData - 1);
\r
768 int destinationLength = 0;
\r
774 destinationLength = 1;
\r
780 destinationLength = 2;
\r
786 destinationLength = 4;
\r
789 // unsupported type for integer destination (float or var-length strings)
\r
793 printf("ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
796 // unknown tag type
\r
798 printf("ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
802 // store in destination
\r
804 memcpy(&destination, pTagData, destinationLength);
\r
808 // tag not found, return failure
\r
813 bool BamAlignment::GetTag(const std::string& tag, int32_t& destination) const {
\r
814 return GetTag(tag, (uint32_t&)destination);
\r
818 bool BamAlignment::GetTag(const std::string& tag, float& destination) const {
\r
820 // make sure tag data exists
\r
821 if ( SupportData.HasCoreOnly || TagData.empty() )
\r
824 // localize the tag data
\r
825 char* pTagData = (char*)TagData.data();
\r
826 const unsigned int tagDataLength = TagData.size();
\r
827 unsigned int numBytesParsed = 0;
\r
829 // if tag found, determine data byte-length, store data in readGroup, return success
\r
830 if ( FindTag(tag, pTagData, tagDataLength, numBytesParsed) ) {
\r
831 //pTagData += numBytesParsed;
\r
833 // determine data byte-length
\r
834 const char type = *(pTagData - 1);
\r
835 int destinationLength = 0;
\r
842 destinationLength = 1;
\r
848 destinationLength = 2;
\r
855 destinationLength = 4;
\r
858 // unsupported type (var-length strings)
\r
861 printf("ERROR: Cannot store tag of type %c in integer destination\n", type);
\r
864 // unknown tag type
\r
866 printf("ERROR: Unknown tag storage class encountered: [%c]\n", type);
\r
870 // store in destination
\r
872 memcpy(&destination, pTagData, destinationLength);
\r
876 // tag not found, return failure
\r
881 bool BamAlignment::RemoveTag(const std::string& tag) {
\r
883 // BamAlignments fetched using BamReader::GetNextAlignmentCore() are not allowed
\r
884 // also, return false if no data present to remove
\r
885 if ( SupportData.HasCoreOnly || TagData.empty() ) return false;
\r
887 // localize the tag data
\r
888 char* pOriginalTagData = (char*)TagData.data();
\r
889 char* pTagData = pOriginalTagData;
\r
890 const unsigned int originalTagDataLength = TagData.size();
\r
891 unsigned int newTagDataLength = 0;
\r
892 unsigned int numBytesParsed = 0;
\r
894 // if tag found, store data in readGroup, return success
\r
895 if ( FindTag(tag, pTagData, originalTagDataLength, numBytesParsed) ) {
\r
897 char newTagData[originalTagDataLength];
\r
899 // copy original tag data up til desired tag
\r
901 numBytesParsed -= 3;
\r
902 const unsigned int beginningTagDataLength = numBytesParsed;
\r
903 newTagDataLength += beginningTagDataLength;
\r
904 memcpy(newTagData, pOriginalTagData, numBytesParsed);
\r
906 // skip to next tag (if tag for removal is last, return true)
\r
907 const char* pTagStorageType = pTagData + 2;
\r
909 numBytesParsed += 3;
\r
910 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return true;
\r
912 // copy everything from current tag (the next one after tag for removal) to end
\r
913 const unsigned int skippedDataLength = (numBytesParsed - beginningTagDataLength);
\r
914 const unsigned int endTagDataLength = originalTagDataLength - beginningTagDataLength - skippedDataLength;
\r
915 memcpy(newTagData + beginningTagDataLength, pTagData, endTagDataLength );
\r
917 // save new tag data
\r
918 TagData.assign(newTagData, beginningTagDataLength + endTagDataLength);
\r
922 // tag not found, no removal - return failure
\r
927 bool BamAlignment::FindTag(const std::string& tag, char* &pTagData, const unsigned int& tagDataLength, unsigned int& numBytesParsed) {
\r
929 while ( numBytesParsed < tagDataLength ) {
\r
931 const char* pTagType = pTagData;
\r
932 const char* pTagStorageType = pTagData + 2;
\r
934 numBytesParsed += 3;
\r
936 // check the current tag, return true on match
\r
937 if ( std::strncmp(pTagType, tag.c_str(), 2) == 0 )
\r
940 // get the storage class and find the next tag
\r
941 if ( *pTagStorageType == '\0' ) return false;
\r
942 if ( !SkipToNextTag(*pTagStorageType, pTagData, numBytesParsed) ) return false;
\r
943 if ( *pTagData == '\0' ) return false;
\r
946 // checked all tags, none match
\r
951 bool BamAlignment::SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {
\r
953 switch(storageType) {
\r
964 numBytesParsed += 2;
\r
971 numBytesParsed += 4;
\r
981 // increment for null-terminator
\r
988 printf("ERROR: Unknown tag storage class encountered: [%c]\n", storageType);
\r
996 } // namespace BamTools
\r