BamAux.h

   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: 31 March 2010 (DB)\r
   7 // ---------------------------------------------------------------------------\r
   8 // Provides the basic constants, data structures, etc. for using BAM files\r
   9 // ***************************************************************************\r
  10 \r
  11 #ifndef BAMAUX_H\r
  12 #define BAMAUX_H\r
  13 \r
  14 // C inclues\r
  15 #include <cstdio>\r
  16 #include <cstdlib>\r
  17 #include <cstring>\r
  18 \r
  19 // C++ includes\r
  20 #include <exception>\r
  21 #include <map>\r
  22 #include <string>\r
  23 #include <utility>\r
  24 #include <vector>\r
  25 \r
  26 // Platform-specific type definitions\r
  27 #ifndef BAMTOOLS_TYPES\r
  28 #define BAMTOOLS_TYPES\r
  29     #ifdef _MSC_VER\r
  30         typedef char                 int8_t;\r
  31         typedef unsigned char       uint8_t;\r
  32         typedef short               int16_t;\r
  33         typedef unsigned short     uint16_t;\r
  34         typedef int                 int32_t;\r
  35         typedef unsigned int       uint32_t;\r
  36         typedef long long           int64_t;\r
  37         typedef unsigned long long uint64_t;\r
  38     #else\r
  39         #include <stdint.h>\r
  40     #endif\r
  41 #endif // BAMTOOLS_TYPES\r
  42 \r
  43 namespace BamTools {\r
  44 \r
  45 // BAM constants\r
  46 const int BAM_CORE_SIZE   = 32;\r
  47 const int BAM_CMATCH      = 0;\r
  48 const int BAM_CINS        = 1;\r
  49 const int BAM_CDEL        = 2;\r
  50 const int BAM_CREF_SKIP   = 3;\r
  51 const int BAM_CSOFT_CLIP  = 4;\r
  52 const int BAM_CHARD_CLIP  = 5;\r
  53 const int BAM_CPAD        = 6;\r
  54 const int BAM_CIGAR_SHIFT = 4;\r
  55 const int BAM_CIGAR_MASK  = ((1 << BAM_CIGAR_SHIFT) - 1);\r
  56 \r
  57 // BAM index constants\r
  58 const int MAX_BIN           = 37450;    // =(8^6-1)/7+1\r
  59 const int BAM_MIN_CHUNK_GAP = 32768;\r
  60 const int BAM_LIDX_SHIFT    = 14;\r
  61 \r
  62 // Explicit variable sizes\r
  63 const int BT_SIZEOF_INT = 4;\r
  64 \r
  65 struct CigarOp;\r
  66 \r
  67 struct BamAlignment {\r
  68 \r
  69     // Queries against alignment flag\r
  70     public:\r
  71         // Returns true if this read is a PCR duplicate (determined by external app)\r
  72         bool IsDuplicate(void) const { return ( (AlignmentFlag & DUPLICATE) != 0 ); }\r
  73         // Returns true if this read failed quality control (determined by external app)\r
  74         bool IsFailedQC(void) const { return ( (AlignmentFlag & QC_FAILED) != 0 ); }\r
  75         // Returns true if alignment is first mate on read\r
  76         bool IsFirstMate(void) const { return ( (AlignmentFlag & READ_1) != 0 ); }\r
  77         // Returns true if alignment is mapped\r
  78         bool IsMapped(void) const { return ( (AlignmentFlag & UNMAPPED) == 0 ); }\r
  79         // Returns true if alignment's mate is mapped\r
  80         bool IsMateMapped(void) const { return ( (AlignmentFlag & MATE_UNMAPPED) == 0 ); }\r
  81         // Returns true if alignment's mate mapped to reverse strand\r
  82         bool IsMateReverseStrand(void) const { return ( (AlignmentFlag & MATE_REVERSE)  != 0 ); }\r
  83         // Returns true if alignment part of paired-end read\r
  84         bool IsPaired(void) const { return ( (AlignmentFlag & PAIRED) != 0 ); }\r
  85         // Returns true if this position is primary alignment (determined by external app)\r
  86         bool IsPrimaryAlignment(void) const  { return ( (AlignmentFlag & SECONDARY) == 0 ); }\r
  87         // Returns true if alignment is part of read that satisfied paired-end resolution (determined by external app)\r
  88         bool IsProperPair(void) const { return ( (AlignmentFlag & PROPER_PAIR) != 0 ); }\r
  89         // Returns true if alignment mapped to reverse strand\r
  90         bool IsReverseStrand(void) const { return ( (AlignmentFlag & REVERSE) != 0 ); }\r
  91         // Returns true if alignment is second mate on read\r
  92         bool IsSecondMate(void) const { return ( (AlignmentFlag & READ_2) != 0 ); }\r
  93 \r
  94     // Manipulate alignment flag\r
  95     public:\r
  96         // Sets "PCR duplicate" bit \r
  97         void SetIsDuplicate(bool ok) { if (ok) AlignmentFlag |= DUPLICATE; else AlignmentFlag &= ~DUPLICATE; }\r
  98         // Sets "failed quality control" bit\r
  99         void SetIsFailedQC(bool ok) { if (ok) AlignmentFlag |= QC_FAILED; else AlignmentFlag &= ~QC_FAILED; }\r
 100         // Sets "alignment is first mate" bit\r
 101         void SetIsFirstMate(bool ok) { if (ok) AlignmentFlag |= READ_1; else AlignmentFlag &= ~READ_1; }\r
 102         // Sets "alignment's mate is mapped" bit\r
 103         void SetIsMateUnmapped(bool ok) { if (ok) AlignmentFlag |= MATE_UNMAPPED; else AlignmentFlag &= ~MATE_UNMAPPED; }\r
 104         // Sets "alignment's mate mapped to reverse strand" bit\r
 105         void SetIsMateReverseStrand(bool ok) { if (ok) AlignmentFlag |= MATE_REVERSE; else AlignmentFlag &= ~MATE_REVERSE; }\r
 106         // Sets "alignment part of paired-end read" bit\r
 107         void SetIsPaired(bool ok) { if (ok) AlignmentFlag |= PAIRED; else AlignmentFlag &= ~PAIRED; }\r
 108         // Sets "alignment is part of read that satisfied paired-end resolution" bit\r
 109         void SetIsProperPair(bool ok) { if (ok) AlignmentFlag |= PROPER_PAIR; else AlignmentFlag &= ~PROPER_PAIR; }\r
 110         // Sets "alignment mapped to reverse strand" bit\r
 111         void SetIsReverseStrand(bool ok) { if (ok) AlignmentFlag |= REVERSE; else AlignmentFlag &= ~REVERSE; }\r
 112         // Sets "position is primary alignment (determined by external app)"\r
 113         void SetIsSecondaryAlignment(bool ok)  { if (ok) AlignmentFlag |= SECONDARY; else AlignmentFlag &= ~SECONDARY; }\r
 114         // Sets "alignment is second mate on read" bit\r
 115         void SetIsSecondMate(bool ok) { if (ok) AlignmentFlag |= READ_2; else AlignmentFlag &= ~READ_2; }\r
 116         // Sets "alignment is mapped" bit\r
 117         void SetIsUnmapped(bool ok) { if (ok) AlignmentFlag |= UNMAPPED; else AlignmentFlag &= ~UNMAPPED; }\r
 118 \r
 119     public:\r
 120 \r
 121         // get "RG" tag data\r
 122         bool GetReadGroup(std::string& readGroup) const {\r
 123 \r
 124             if ( TagData.empty() ) { return false; }\r
 125 \r
 126             // localize the tag data\r
 127             char* pTagData = (char*)TagData.data();\r
 128             const unsigned int tagDataLen = TagData.size();\r
 129             unsigned int numBytesParsed = 0;\r
 130 \r
 131             bool foundReadGroupTag = false;\r
 132             while( numBytesParsed < tagDataLen ) {\r
 133 \r
 134                 const char* pTagType = pTagData;\r
 135                 const char* pTagStorageType = pTagData + 2;\r
 136                 pTagData       += 3;\r
 137                 numBytesParsed += 3;\r
 138 \r
 139                 // check the current tag\r
 140                 if ( std::strncmp(pTagType, "RG", 2) == 0 ) {\r
 141                     foundReadGroupTag = true;\r
 142                     break;\r
 143                 }\r
 144 \r
 145                 // get the storage class and find the next tag\r
 146                 if (*pTagStorageType == '\0') { return false; }\r
 147                 SkipToNextTag( *pTagStorageType, pTagData, numBytesParsed );\r
 148                 if (*pTagData == '\0') { return false; }\r
 149             }\r
 150 \r
 151             // return if the read group tag was not present\r
 152             if ( !foundReadGroupTag ) { return false; }\r
 153 \r
 154             // assign the read group\r
 155             const unsigned int readGroupLen = std::strlen(pTagData);\r
 156             readGroup.resize(readGroupLen);\r
 157             std::memcpy( (char*)readGroup.data(), pTagData, readGroupLen );\r
 158             return true;\r
 159         }\r
 160 \r
 161         // get "NM" tag data - contributed by Aaron Quinlan\r
 162         bool GetEditDistance(uint8_t& editDistance) const {\r
 163 \r
 164             if ( TagData.empty() ) { return false; }\r
 165 \r
 166             // localize the tag data\r
 167             char* pTagData = (char*)TagData.data();\r
 168             const unsigned int tagDataLen = TagData.size();\r
 169             unsigned int numBytesParsed = 0;\r
 170 \r
 171             bool foundEditDistanceTag = false;\r
 172             while( numBytesParsed < tagDataLen ) {\r
 173 \r
 174                 const char* pTagType = pTagData;\r
 175                 const char* pTagStorageType = pTagData + 2;\r
 176                 pTagData       += 3;\r
 177                 numBytesParsed += 3;\r
 178 \r
 179                 // check the current tag\r
 180                 if ( strncmp(pTagType, "NM", 2) == 0 ) {\r
 181                     foundEditDistanceTag = true;\r
 182                     break;\r
 183                 }\r
 184 \r
 185                 // get the storage class and find the next tag\r
 186                 if (*pTagStorageType == '\0') { return false; }\r
 187                 SkipToNextTag( *pTagStorageType, pTagData, numBytesParsed );\r
 188                 if (*pTagData == '\0') { return false; }\r
 189             }\r
 190             // return if the edit distance tag was not present\r
 191             if ( !foundEditDistanceTag ) { return false; }\r
 192 \r
 193             // assign the editDistance value\r
 194             std::memcpy(&editDistance, pTagData, 1);\r
 195             return true;\r
 196         }\r
 197 \r
 198     private:\r
 199         static void SkipToNextTag(const char storageType, char* &pTagData, unsigned int& numBytesParsed) {\r
 200             switch(storageType) {\r
 201 \r
 202                 case 'A':\r
 203                 case 'c':\r
 204                 case 'C':\r
 205                         ++numBytesParsed;\r
 206                         ++pTagData;\r
 207                         break;\r
 208 \r
 209                 case 's':\r
 210                 case 'S':\r
 211                         numBytesParsed += 2;\r
 212                         pTagData       += 2;\r
 213                         break;\r
 214 \r
 215                 case 'f':\r
 216                 case 'i':\r
 217                 case 'I':\r
 218                         numBytesParsed += 4;\r
 219                         pTagData       += 4;\r
 220                         break;\r
 221 \r
 222                 case 'Z':\r
 223                 case 'H':\r
 224                         while(*pTagData) {\r
 225                             ++numBytesParsed;\r
 226                             ++pTagData;\r
 227                         }\r
 228                         // ---------------------------\r
 229                         // Added: 3-25-2010 DWB\r
 230                         // Contributed: ARQ\r
 231                         // Fixed: error parsing variable length tag data\r
 232                         ++pTagData;\r
 233                         // ---------------------------\r
 234                         break;\r
 235 \r
 236                 default:\r
 237                         printf("ERROR: Unknown tag storage class encountered: [%c]\n", *pTagData);\r
 238                         exit(1);\r
 239             }\r
 240         }\r
 241 \r
 242     // Data members\r
 243     public:\r
 244         std::string  Name;              // Read name\r
 245         int32_t      Length;            // Query length\r
 246         std::string  QueryBases;        // 'Original' sequence (as reported from sequencing machine)\r
 247         std::string  AlignedBases;      // 'Aligned' sequence (includes any indels, padding, clipping)\r
 248         std::string  Qualities;         // FASTQ qualities (ASCII characters, not numeric values)\r
 249         std::string  TagData;           // Tag data (accessor methods will pull the requested information out)\r
 250         int32_t      RefID;             // ID number for reference sequence\r
 251         int32_t      Position;          // Position (0-based) where alignment starts\r
 252         uint16_t     Bin;               // Bin in BAM file where this alignment resides\r
 253         uint16_t     MapQuality;        // Mapping quality score\r
 254         uint32_t     AlignmentFlag;     // Alignment bit-flag - see Is<something>() methods to query this value, SetIs<something>() methods to manipulate \r
 255         std::vector<CigarOp> CigarData; // CIGAR operations for this alignment\r
 256         int32_t      MateRefID;         // ID number for reference sequence where alignment's mate was aligned\r
 257         int32_t      MatePosition;      // Position (0-based) where alignment's mate starts\r
 258         int32_t      InsertSize;        // Mate-pair insert size\r
 259 \r
 260     // Alignment flag query constants\r
 261     private:\r
 262         enum { PAIRED        = 1,\r
 263                PROPER_PAIR   = 2,\r
 264                UNMAPPED      = 4,\r
 265                MATE_UNMAPPED = 8,\r
 266                REVERSE       = 16,\r
 267                MATE_REVERSE  = 32,\r
 268                READ_1        = 64,\r
 269                READ_2        = 128,\r
 270                SECONDARY     = 256,\r
 271                QC_FAILED     = 512,\r
 272                DUPLICATE     = 1024\r
 273              };\r
 274 };\r
 275 \r
 276 // ----------------------------------------------------------------\r
 277 // Auxiliary data structs & typedefs\r
 278 \r
 279 struct CigarOp {\r
 280     char     Type;   // Operation type (MIDNSHP)\r
 281     uint32_t Length; // Operation length (number of bases)\r
 282 };\r
 283 \r
 284 struct RefData {\r
 285     // data members\r
 286     std::string RefName;          // Name of reference sequence\r
 287     int32_t     RefLength;        // Length of reference sequence\r
 288     bool        RefHasAlignments; // True if BAM file contains alignments mapped to reference sequence\r
 289     // constructor\r
 290     RefData(void)\r
 291         : RefLength(0)\r
 292         , RefHasAlignments(false)\r
 293     { }\r
 294 };\r
 295 \r
 296 typedef std::vector<RefData>      RefVector;\r
 297 typedef std::vector<BamAlignment> BamAlignmentVector;\r
 298 \r
 299 // ----------------------------------------------------------------\r
 300 // Indexing structs & typedefs\r
 301 \r
 302 struct Chunk {\r
 303     // data members\r
 304     uint64_t Start;\r
 305     uint64_t Stop;\r
 306     // constructor\r
 307     Chunk(const uint64_t& start = 0, const uint64_t& stop = 0)\r
 308         : Start(start)\r
 309         , Stop(stop)\r
 310     { }\r
 311 };\r
 312 \r
 313 inline\r
 314 bool ChunkLessThan(const Chunk& lhs, const Chunk& rhs) {\r
 315     return lhs.Start < rhs.Start;\r
 316 }\r
 317 \r
 318 typedef std::vector<Chunk> ChunkVector;\r
 319 typedef std::map<uint32_t, ChunkVector> BamBinMap;\r
 320 typedef std::vector<uint64_t> LinearOffsetVector;\r
 321 \r
 322 struct ReferenceIndex {\r
 323     // data members\r
 324     BamBinMap Bins;\r
 325     LinearOffsetVector Offsets;\r
 326     // constructor\r
 327     ReferenceIndex(const BamBinMap& binMap = BamBinMap(),\r
 328                    const LinearOffsetVector& offsets = LinearOffsetVector())\r
 329         : Bins(binMap)\r
 330         , Offsets(offsets)\r
 331     { }\r
 332 };\r
 333 \r
 334 typedef std::vector<ReferenceIndex> BamIndex;\r
 335 \r
 336 // ----------------------------------------------------------------\r
 337 // Added: 3-35-2010 DWB\r
 338 // Fixed: Routines to provide endian-correctness\r
 339 // ----------------------------------------------------------------\r
 340 \r
 341 // returns true if system is big endian\r
 342 inline bool SystemIsBigEndian(void) {\r
 343    const uint16_t one = 0x0001;\r
 344    return ((*(char*) &one) == 0 );\r
 345 }\r
 346 \r
 347 // swaps endianness of 16-bit value 'in place'\r
 348 inline void SwapEndian_16(int16_t& x) {\r
 349     x = ((x >> 8) | (x << 8));\r
 350 }\r
 351 \r
 352 inline void SwapEndian_16(uint16_t& x) {\r
 353     x = ((x >> 8) | (x << 8));\r
 354 }\r
 355 \r
 356 // swaps endianness of 32-bit value 'in-place'\r
 357 inline void SwapEndian_32(int32_t& x) {\r
 358     x = ( (x >> 24) | \r
 359          ((x << 8) & 0x00FF0000) | \r
 360          ((x >> 8) & 0x0000FF00) | \r
 361           (x << 24)\r
 362         );\r
 363 }\r
 364 \r
 365 inline void SwapEndian_32(uint32_t& x) {\r
 366     x = ( (x >> 24) | \r
 367          ((x << 8) & 0x00FF0000) | \r
 368          ((x >> 8) & 0x0000FF00) | \r
 369           (x << 24)\r
 370         );\r
 371 }\r
 372 \r
 373 // swaps endianness of 64-bit value 'in-place'\r
 374 inline void SwapEndian_64(int64_t& x) {\r
 375     x = ( (x >> 56) | \r
 376          ((x << 40) & 0x00FF000000000000ll) |\r
 377          ((x << 24) & 0x0000FF0000000000ll) |\r
 378          ((x << 8)  & 0x000000FF00000000ll) |\r
 379          ((x >> 8)  & 0x00000000FF000000ll) |\r
 380          ((x >> 24) & 0x0000000000FF0000ll) |\r
 381          ((x >> 40) & 0x000000000000FF00ll) |\r
 382           (x << 56)\r
 383         );\r
 384 }\r
 385 \r
 386 inline void SwapEndian_64(uint64_t& x) {\r
 387     x = ( (x >> 56) | \r
 388          ((x << 40) & 0x00FF000000000000ll) |\r
 389          ((x << 24) & 0x0000FF0000000000ll) |\r
 390          ((x << 8)  & 0x000000FF00000000ll) |\r
 391          ((x >> 8)  & 0x00000000FF000000ll) |\r
 392          ((x >> 24) & 0x0000000000FF0000ll) |\r
 393          ((x >> 40) & 0x000000000000FF00ll) |\r
 394           (x << 56)\r
 395         );\r
 396 }\r
 397 \r
 398 inline void SwapEndian_16p(char* data) {\r
 399     uint16_t& value = (uint16_t&)*data; \r
 400     SwapEndian_16(value);\r
 401 }\r
 402 \r
 403 inline void SwapEndian_32p(char* data) {\r
 404     uint32_t& value = (uint32_t&)*data; \r
 405     SwapEndian_32(value);\r
 406 }\r
 407 \r
 408 inline void SwapEndian_64p(char* data) {\r
 409     uint64_t& value = (uint64_t&)*data; \r
 410     SwapEndian_64(value);\r
 411 }\r
 412 \r
 413 } // namespace BamTools\r
 414 \r
 415 #endif // BAMAUX_H\r